This study is a monocentric, prospective, observational investigation conducted in Nancy Hospital between 10.10.2017 and 24.05.2022, enrolling 50 adult patients with refractory cardiogenic shock requiring veno-arterial extracorporeal membrane oxygenation (VA-ECMO) support (NCT03327493).
The study was carried out in accordance with the Declaration of Helsinki and approved by the local institutional ethics committee. Written informed consent was obtained from all patients or their legal representatives prior to inclusion.
Cardiogenic shock was defined as a systolic blood pressure (SBP) <90 mmHg or a mean arterial pressure (MAP) <65 mmHg despite adequate volume resuscitation, accompanied by signs of peripheral hypoperfusion (e.g., cold extremities, oliguria, altered mental status) and a cardiac index <2.2 L/min/m². Refractoriness was defined as a hyporesponsiveness to norepinephrine and/or the persistence of profound clinical signs of hypoperfusion despite optimal resuscitation efforts. Eligible patients were adults (≥18 years old) affiliated with a health insurance system. Cardiogenic shock etiologies included acute ischemic coronary artery disease, primary or ischemic-related dilated cardiomyopathy, myocarditis and stress-induced cardiomyopathy. Patients were excluded if cardiogenic shock was secondary to cardiotoxic poisoning, if they were pregnant, or were under legal protective supervision (e.g., court-appointed guardianship).
Clinical, laboratory and echocardiographic parameters were prospectively collected at three predefined timepoints:Day 0 (ECMO implantation), Day 3–5 (early course) and Day of ECMO weaning (ECMO explantation). Blood samples were drawn at each timepoint, immediately centrifuged and plasma was stored at –80°C until analysis. Samples were then transferred to the XXX for cells analysis and INSERM UMR942 laboratory (Paris, France) for cytokine analysis. Quantification of adrenoreceptors β1 and β2 on monocytes and CD4⁺ T lymphocytes was performed using flow cytometry. Cytokine concentrations in plasma, including interleukin (IL) -33, IL-1β, IL-2, IL-4, IL-5, IL-6, IL-8/CXCL8, IL-10, GDF-15, IL-12p70, IL-17A/CTLA-8, IFN-γ, TNF-α and granzyme B, were measured using [XXX].
The expression of β1 (ADRB1) and β2 (ADRB2) adrenergic receptors on circulating immune cell, namely CD4+ T lymphocytes and monocytes, was assessed via flow cytometry. XXX [to be completed from Manon]. Results were expressed as the percentage of cells that were positive for the respective expression receptors.
##Cytokine Quantification
The aforementioned panel of circulating cytokines and immune mediators was quantified from plasma samples. These biomarkers were selected for their roles in pro- and anti-inflammatory signaling. Measurements were performed using a multiplex bead-based immunoassay (Luminex platform), following the manufacturer’s instructions. XXX [to be completed from Malha].
##ECMO Weaning Protocol
A standardized weaning protocol was applied to all patients prior to ECMO explantation. At a constant mean arterial pressure (MAP) between 65 and 75 mmHg, cardiac output was assessed by transthoracic echocardiography at two predefined timepoints: 1. Baseline under full support (ECLS flow at 3.0 L/min) and 2. Low-flow test after 45 minutes of reduced ECMO flow (1.5 L/min). Measurements included velocity time integral (VTI) at the left ventricular outflow tract and left ventricular ejection fraction. Changes in norepinephrine and dobutamin dosage during this period were also documented. The decision to explant ECMO was based on clinical tolerance, hemodynamic stability and echocardiographic parameters during the weaning trial.
#Data load
# Blood sample at day 0 before ECMO implant?
df %>%
select(J0_date, j0_bilan_date) %>%
distinct() %>%
slice(1:10)
# A tibble: 10 × 2
J0_date j0_bilan_date
<chr> <chr>
1 16.10.2017 17.10.2017
2 25.10.2017 25.10.2017
3 17.01.2018 17.01.2018
4 25.01.2018 25.01.2018
5 06.02.2018 06.02.2018
6 20.02.2018 20.02.2018
7 02.03.2018 02.03.2018
8 30.05.2018 30.05.2018
9 14.06.2018 14.06.2018
10 27.06.2018 27.06.2018
df <- df %>%
mutate(
J0_date = as.Date(J0_date, format = "%d.%m.%Y"),
j0_bilan_date = as.Date(j0_bilan_date, format = "%d.%m.%Y"),
bilan_after_ecmo = case_when(
is.na(J0_date) | is.na(j0_bilan_date) ~ NA_character_,
j0_bilan_date > J0_date ~ "Yes",
TRUE ~ "No"
)
)
table(df$bilan_after_ecmo, useNA = "ifany")
No Yes
48 1
## 2 patients had blood sample AFTER ECMO implant.
df %>%
filter(bilan_after_ecmo == "Yes") %>%
select(ID, j0_bilan_date, J0_date, bilan_after_ecmo)
# A tibble: 1 × 4
ID j0_bilan_date J0_date bilan_after_ecmo
<chr> <date> <date> <chr>
1 1-FG 2017-10-17 2017-10-16 Yes
# Patient 1 (2017-10-17 > 2017-10-16) and 15 (2019-04-19 > 2018-04-19). Patients 15 has probably an error. In PAtinet 1 the blood sample was after 12 hours.
df <- df %>%
# Applica solo alle righe dove bilan_after_ecmo è "No"
mutate(
heure_check = case_when(
bilan_after_ecmo == "No" &
!is.na(J0_heure) & !is.na(j0_bilan_heure) ~ if_else(
j0_bilan_heure > J0_heure, "Bilan after ECMO (same day)",
"Bilan before/same as ECMO (same day)"
),
TRUE ~ NA_character_
)
)
table(df$bilan_after_ecmo, useNA = "ifany")
No Yes
48 1
df %>%
filter(heure_check == "Bilan after ECMO (same day)") %>%
select(ID, J0_date, J0_heure, j0_bilan_date, j0_bilan_heure, heure_check)
# A tibble: 0 × 6
# ℹ 6 variables: ID <chr>, J0_date <date>, J0_heure <time>,
# j0_bilan_date <date>, j0_bilan_heure <time>, heure_check <chr>
# Among patients with laboratory tests on the same day of ECMO initiation, none had blood sampling performed after ECMO implantation.
# Mean time between blood sample and ECMO implant, patients 1 e 15 excluded.
df_clean <- df %>%
filter(!str_starts(ID, "1-") & !str_starts(ID, "15-")) %>%
mutate(
j0_bilan_date = as.Date(j0_bilan_date, format = "%d.%m.%Y"),
J0_date = as.Date(J0_date, format = "%d.%m.%Y"),
datetime_bilan = as.POSIXct(paste(j0_bilan_date, as.character(j0_bilan_heure)), format = "%Y-%m-%d %H:%M", tz = "UTC"),
datetime_ecmo = as.POSIXct(paste(J0_date, as.character(J0_heure)), format = "%Y-%m-%d %H:%M", tz = "UTC"),
time_diff_mins = as.numeric(difftime(datetime_bilan, datetime_ecmo, units = "mins"))
) %>%
filter(!is.na(time_diff_mins))
df_clean %>%
summarise(
mean_minutes = mean(time_diff_mins),
sd_minutes = sd(time_diff_mins),
min = min(time_diff_mins),
max = max(time_diff_mins)
)
# A tibble: 1 × 4
mean_minutes sd_minutes min max
<dbl> <dbl> <dbl> <dbl>
1 -205. 145. -660 0
In 2 patients, the day 0 laboratory test occurred after ECMO initiation. Day-0 blood samples were collected on average 3.4 hours (SD 2.4 hours) before ECMO implantation, ranging from 11 hours before to the exact time of ECMO cannulation.
df <- df %>%
mutate(
# Basic formatting
Temp_D0 = round(Temp_D0, 1),
Age = as.numeric(Age),
Gender = case_when(
Gender == 0 ~ "Men",
Gender == 1 ~ "Women",
TRUE ~ as.character(Gender)
),
Preecmo_tte_diam_LV = Preecmo_tte_diam_LV * 10,
VIS_D0 = (DOBUcum_D0 * 1000) / (Weight * 1440) +
100 * (NADcum_D0 * 1000) / (Weight * 1440),
# Outcome recoding
Outcome = case_when(
Outcome == 1 ~ "ECMO Weaning",
Outcome == 2 ~ "Bridge to Transplant",
Outcome == 3 ~ "Bridge to LVAD",
Outcome == 4 ~ "Death",
TRUE ~ as.character(Outcome)
),
Pulsepressure_D0 = PAS_D0 - PAD_D0,
Outcome_censored = as.numeric(case_when(
Outcome == "Death" ~ "1",
Outcome %in% c("ECMO Weaning", "Bridge to Transplant", "Bridge to LVAD") ~ "0"
)),
outcome_date = str_replace_all(outcome_date, "\\.", "/"),
icu_admiss_date = str_replace_all(icu_admiss_date, "\\.", "/"),
icu_discharge_date = str_replace_all(icu_discharge_date, "\\.", "/"),
outcome_date = dmy(outcome_date),
icu_admiss_date = dmy(icu_admiss_date),
icu_discharge_date = dmy(icu_discharge_date),
date_levo = dmy(date_levo),
ecmo_start_date = dmy(ecmo_start_date),
ecmo_stop_date = dmy(ecmo_stop_date),
diff_days = time_length(interval(icu_admiss_date, outcome_date), "days"),
Time_hosp = time_length(interval(icu_admiss_date, icu_discharge_date), "days"),
as.numeric(icu_discharge_date - icu_admiss_date),
Time_ecmo = time_length(interval(ecmo_start_date, ecmo_stop_date), "days"),
Time_levo = time_length(interval(ecmo_start_date, date_levo), "days"),
diff_days_J28 = case_when(
diff_days > 28 ~ 28,
TRUE ~ diff_days
),
Outcome_censored_28 = case_when(
diff_days > 28 & Outcome_censored == 1 ~ 0,
TRUE ~ Outcome_censored
),
ADR1_D0_median = case_when(
ADR1_D0 > median(ADR1_D0, na.rm = TRUE) ~ "High",
ADR1_D0 <= median(ADR1_D0, na.rm = TRUE) ~ "Low"
),
ADR2_D0_median = case_when(
ADR2_D0 > median(ADR2_D0, na.rm = TRUE) ~ "High",
ADR2_D0 <= median(ADR2_D0, na.rm = TRUE) ~ "Low"
),
ADR1m_D0_median = case_when(
ADR1m_J0 > median(ADR1m_J0, na.rm = TRUE) ~ "High",
ADR1m_J0 <= median(ADR1m_J0, na.rm = TRUE) ~ "Low",
TRUE ~ NA_character_
),
ADR2m_D0_median = case_when(
ADR2m_J0 > median(ADR2m_J0, na.rm = TRUE) ~ "High",
ADR2m_J0 <= median(ADR2m_J0, na.rm = TRUE) ~ "Low",
TRUE ~ NA_character_
),
Betablocker = case_when(
Betablocker == 1 ~ "YES",
Betablocker == 0 ~ "NO"
),
dose_tot_dobu_yes_no = droplevels(as.factor(case_when(
dose_tot_dobu > 0 ~ "YES",
dose_tot_dobu == 0 ~ "NO"
))),
levo_in_icu = case_when(
Levosimendan == 1 & Time_levo >= 0 ~ "Yes",
Levosimendan == 1 & Time_levo < 0 ~ "No",
Levosimendan == 0 ~ "No"
),
Levosimendan = case_when(
Levosimendan == 1 ~ "YES",
Levosimendan == 0 ~ "NO",
TRUE ~ NA_character_
)
)
var.cat <- c("Gender", "HTN", "DM", "Immunodepression", "CKD", "Neurologic_deficit", "Chronic_respiratory_disease", "HF", "Cardiac_arrest_before_canul", "IABP_D0", "Intubation", "EER_D0", "Outcome", "Betablocker", "Alphablocker", "Levosimendan", "Death")
var.cont <- c("Age", "BMI", "PAS_D0", "PAD_D0", "PAM_D0", "Pulsepressure_D0", "HR_D0", "Temp_D0", "Preecmo_tte_ef", "Preecmo_tte_vtiao", "Preecmo_tte_diam_LV", "Preecmo_tte_tapse", "Preecmo_tte_vmax_s_mit_lat_LV", "NADcum_D0", "DOBUcum_D0", "VIS_D0", "pH_D0", "Sao2_D0", "Paco2_D0", "Lact_D0", "Hco3_D0", "Creat_D0", "Urea_D0", "ALAT_D0", "ASAT_D0", "Bili_D0", "PT_D0", "Tropo_i_hs_D0", "Ntprobnp_D0", "DPP3_D0", "Plq_D0", "Lenght_eer", "Lenght_vm", "Time_hosp", "Time_ecmo")
var.tot <- c(var.cat, var.cont)
table1
Variables | N | Overall |
|---|---|---|
Age (years) | 50 | 59.00 [49.00, 66.00] |
Gender = Female (%) | 50 | 8 (16.3) |
Body Mass Index (kg/m²) | 47 | 26.98 [23.72, 30.00] |
Hypertension (%) | 50 | 19 (38.8) |
Diabetes Mellitus (%) | 50 | 10 (20.4) |
Heart Failure (%) | 50 | 17 (34.7) |
Chronic Kidney Disease (%) | 50 | 5 (10.2) |
Chronic Respiratory Disease (%) | 50 | 3 ( 6.1) |
Neurological Deficit (%) | 50 | 4 ( 8.2) |
Immunosuppression (%) | 50 | 4 ( 8.2) |
Systolic Blood Pressure (mmHg) | 50 | 96.00 [86.00, 109.00] |
Diastolic Blood Pressure (mmHg) | 50 | 63.00 [56.00, 72.00] |
Mean Arterial Pressure (mmHg) | 50 | 72.00 [68.00, 81.00] |
Pulse Pressure (mmHg) | 50 | 31.00 [22.00, 47.00] |
Heart Rate (bpm) | 50 | 96.00 [86.00, 118.00] |
Temperature (°C) | 50 | 36.00 [35.30, 36.80] |
Cardiac Arrest Before Cannulation (%) | 50 | 27 (55.1) |
Intra-Aortic Balloon Pump (%) | 50 | 23 (46.9) |
Intubation (%) | 50 | 36 (73.5) |
Renal Replacement Therapy (%) | 50 | 7 (14.3) |
Pre-ECMO Ejection Fraction (%) | 46 | 20.00 [15.00, 25.00] |
Pre-ECMO Aortic Velocity Time Integral (cm) | 41 | 7.00 [5.27, 8.25] |
Pre-ECMO Left Ventricular Diameter (mm) | 31 | 59.00 [48.50, 65.25] |
Pre-ECMO TAPSE (mm) | 38 | 12.50 [8.00, 16.00] |
Pre-ECMO Mitral Lateral Velocity (cm/s) | 26 | 4.50 [4.00, 5.80] |
Beta-blocker (%) | 50 | 22 (44.9) |
Alpha-blocker (%) | 50 | 3 ( 6.1) |
Cumulative Norepinephrine Dose (mg) | 50 | 20.00 [0.00, 86.50] |
Cumulative Dobutamine Dose (mg) | 50 | 200.00 [0.00, 510.00] |
Vasoactive-Inotropic Score (VIS) at Day 0 | 48 | 18.52 [4.86, 58.75] |
Levosimendan (%) | 50 | 29 (59.2) |
pH | 50 | 7.42 [7.31, 7.49] |
Oxygen Saturation (%) | 49 | 98.00 [96.00, 99.00] |
PaCO2 (mmHg) | 50 | 32.70 [25.90, 37.60] |
Lactate (mmol/L) | 49 | 2.75 [1.70, 4.70] |
Bicarbonate (mmol/L) | 50 | 19.00 [16.80, 23.00] |
Creatinine (µmol/L) | 50 | 149.00 [104.00, 202.00] |
Urea (mg/dL) | 49 | 11.21 [8.14, 18.86] |
ALT (IU/L) | 50 | 264.00 [91.00, 804.00] |
AST (IU/L) | 50 | 656.00 [171.00, 1202.00] |
Bilirubin (mg/dL) | 49 | 21.00 [12.00, 38.50] |
Prothrombin Time (s) | 38 | 56.00 [37.00, 73.00] |
High Sensitivity Troponin I (ng/L) | 35 | 34363.50 [383.25, 186316.00] |
NT-proBNP (pg/mL) | 24 | 10882.00 [2248.00, 26995.00] |
DPP3 (ng/mL) | 50 | 72.90 [50.00, 164.10] |
Platelet Count (10⁹/L) | 50 | 176.00 [118.00, 271.00] |
Outcome (%) | 50 |
|
Death in ICU | 50 | 22 (44.9) |
Bridge to LVAD | 50 | 4 ( 8.2) |
Bridge to Transplant | 50 | 2 ( 4.1) |
ECMO Weaning | 50 | 21 (42.9) |
Duration of Renal Replacement Therapy (days) | 50 | 0.00 [0.00, 5.00] |
Duration of Mechanical Ventilation (days) | 50 | 5.00 [2.00, 10.00] |
ECMO Duration (days) | 50 | 6.00 [3.00, 9.00] |
ICU Length of Stay (days) | 50 | 14.00 [6.00, 19.00] |
28-day Mortality (%) | 50 | 28 (57.1) |
50 | 49 |
Baseline clinical, demographic and hemodynamic characteristics were summarized for the overall study population (n = 50). One patient was excluded from the analysis due to being under legal guardianship. Categorical variables were reported as counts and percentages, while continuous variables were expressed as medians and interquartile ranges (IQR). The normality of continuous variables was assessed visually using histograms.
for (var in var.cont) {
if (var %in% colnames(T1)) {
print(
ggplot(T1, aes_string(var)) +
geom_histogram(bins = 15, fill = "blue", alpha = 0.5) +
ggtitle(paste("Distribuzione di", var))
)
} else {
message(paste("❌ Variabile", var, "non trovata in T1"))
}
}
table2a
# A tibble: 8 × 2
Parameter `Median [IQR]`
<chr> <chr>
1 ECMO flow (l) 2.4 [1.6–3], n=33
2 Dobutamine cumulative (mg) 0 [0–50], n=33
3 Noradrenaline cumulative (mg) 0 [0–6.2], n=33
4 MAP_baseline (mmHg) 75 [72–80], n=33
5 HR (bpm) 100 [84–103], n=33
6 FEVG (%) 35 [30–43.2], n=32
7 Mitral S' (cm/s) 7.4 [6–9], n=31
8 Lactate (mmol/L) 1.1 [0.9–1.4], n=32
# Create new variables for flow condition and unify variable names
df_long <- df %>%
select(
ID,
jsev_epr_sev_t0_3l_itv,
jsev_epr_sev_t0_3l_pam,
jsev_epr_sev_t0_3l_NAd,
`jsev_epr_sev_t1_1,5l_itv`,
`jsev_epr_sev_t1_1,5l_pam`,
`jsev_epr_sev_t1_1,5l_NAd`,
Weight,
) %>%
mutate(jsev_epr_sev_t0_3l_NAd_w = jsev_epr_sev_t0_3l_NAd/Weight,
`jsev_epr_sev_t1_1,5l_NAd_w` = `jsev_epr_sev_t1_1,5l_NAd`/Weight)%>%
rename(
`VTI (cm)` = jsev_epr_sev_t0_3l_itv,
`MAP (mmHg)` = jsev_epr_sev_t0_3l_pam,
`Noradrenalin (µg/kg/h)` = jsev_epr_sev_t0_3l_NAd_w,
`VTI_1.5L (cm)` = `jsev_epr_sev_t1_1,5l_itv`,
`MAP_1.5L (mmHg)` = `jsev_epr_sev_t1_1,5l_pam`,
`Noradrenalin_1.5L (µg/kg/h)` = `jsev_epr_sev_t1_1,5l_NAd_w`
) %>%
pivot_longer(
cols = -c(ID,jsev_epr_sev_t0_3l_NAd,`jsev_epr_sev_t1_1,5l_NAd`,Weight),
names_to = "Parameter_raw",
values_to = "Value"
) %>%
mutate(
Value = as.numeric(Value),
Condition = case_when(
str_detect(Parameter_raw, "1.5L") ~ "1.5L",
TRUE ~ "3L"
),
Condition = factor(Condition, levels = c("3L", "1.5L")),
Parameter = str_replace(Parameter_raw, "_1.5L", "")
)
# Create summary table
table2b <- df_long %>%
group_by(Parameter, Condition) %>%
summarise(
n = sum(!is.na(Value)),
Median = round(median(Value, na.rm = TRUE), 2),
Q1 = round(quantile(Value, 0.25, na.rm = TRUE), 2),
Q3 = round(quantile(Value, 0.75, na.rm = TRUE), 2),
.groups = "drop"
) %>%
mutate(
Value_IQR = paste0(Median, " [", Q1, "–", Q3, "]; n=", n)
) %>%
select(Parameter, Condition, Value_IQR) %>%
pivot_wider(names_from = Condition, values_from = Value_IQR)
table2b
# A tibble: 3 × 3
Parameter `3L` `1.5L`
<chr> <chr> <chr>
1 MAP (mmHg) 75.5 [70.75–83.5]; n=20 72.5 [70–80.5]; n=28
2 Noradrenalin (µg/kg/h) 0 [0–0]; n=22 0 [0–0]; n=29
3 VTI (cm) 13 [12–15]; n=21 14 [12–16]; n=29
# Wilcoxon test for a parameter
wilcox.test(Value ~ Condition,
data = df_long %>% filter(Parameter == "VTI (cm)"))
Wilcoxon rank sum test with continuity correction
data: Value by Condition
W = 271.5, p-value = 0.5204
alternative hypothesis: true location shift is not equal to 0
No significant difference between VTI at 3L and 1.5L flow conditions at weaning day (Wilcoxon test, p = 0.5).
F1 <- df %>%
select(ID,ADR1_J0, ADR1_J3_J5, ADR1_JS, Outcome)%>%
mutate(Outcomes = factor(case_when(Outcome=="Bridge to LVAD"|Outcome=="Bridge to Transplant"~ "Bridge to transplant or LVAD",T~Outcome), levels=c("Death","Bridge to transplant or LVAD","ECMO Weaning")))%>%
pivot_longer(
cols = c(ADR1_J0, ADR1_J3_J5, ADR1_JS),
names_to = "ADR1_time",
values_to = "value"
)%>%
mutate(time = as.numeric(case_when(ADR1_time=="ADR1_J0"~1,
ADR1_time=="ADR1_J3_J5"~2,
ADR1_time=="ADR1_JS"~3)),
value_log = log10(value))
count_data <- F1 %>%
filter(!is.na(value))%>%
group_by(ADR1_time, Outcomes) %>%
summarise(n = n(), .groups = "drop")%>%
filter(!is.na(Outcomes))
# ===========================
# Modélisation
# ===========================
# Modèle avec interaction
mod1 <- lmerTest::lmer(
value_log ~ time * Outcomes + (1 + time | ID),
data = F1
)
# Résumé du modèle
summary(mod1)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: value_log ~ time * Outcomes + (1 + time | ID)
Data: F1
REML criterion at convergence: 100.8
Scaled residuals:
Min 1Q Median 3Q Max
-2.35927 -0.53586 0.03986 0.52500 1.90217
Random effects:
Groups Name Variance Std.Dev. Corr
ID (Intercept) 0.11537 0.3397
time 0.01505 0.1227 -0.71
Residual 0.09453 0.3075
Number of obs: 93, groups: ID, 46
Fixed effects:
Estimate Std. Error df t value
(Intercept) -0.15452 0.16111 61.94485 -0.959
time 0.12703 0.08955 54.66302 1.419
OutcomesBridge to transplant or LVAD 0.09504 0.28976 38.98471 0.328
OutcomesECMO Weaning 0.32514 0.21202 48.22743 1.534
time:OutcomesBridge to transplant or LVAD 0.09162 0.13978 34.93053 0.655
time:OutcomesECMO Weaning -0.07095 0.10946 43.99886 -0.648
Pr(>|t|)
(Intercept) 0.341
time 0.162
OutcomesBridge to transplant or LVAD 0.745
OutcomesECMO Weaning 0.132
time:OutcomesBridge to transplant or LVAD 0.516
time:OutcomesECMO Weaning 0.520
Correlation of Fixed Effects:
(Intr) time OBttoL OECMOW t:ttoL
time -0.877
OtcBttoLVAD -0.556 0.487
OtcmsECMOWn -0.760 0.666 0.422
t:OBttoLVAD 0.562 -0.641 -0.860 -0.427
tm:OtcECMOW 0.717 -0.818 -0.399 -0.860 0.524
# Modèle sans interaction (modèle réduit)
mod_no_interaction <- lmerTest::lmer(
value_log ~ time + Outcomes + (1 + time | ID),
data = F1
)
# Comparaison des modèles (test de l'interaction)
interact <- anova(mod_no_interaction, mod1)
# Interprétation : p = 0.58 → interaction non significative
# ===========================
# Tests des effets fixes
# ===========================
# ANOVA avec approximation de Kenward-Roger pour les p-values
a <- anova(mod_no_interaction, ddf = "Kenward-Roger")
print(a)
Type III Analysis of Variance Table with Kenward-Roger's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
time 0.48396 0.48396 1 32.800 5.0378 0.03167 *
Outcomes 0.46031 0.23015 2 36.375 2.3938 0.10552
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ===========================
# Vérification des hypothèses
# ===========================
# 1. Diagrammes de diagnostic généraux
plot(mod1) # résidus vs prédictions par groupe
# 2. Résidus vs valeurs ajustées (homoscédasticité)
plot(fitted(mod1), residuals(mod1),
xlab = "Valeurs ajustées", ylab = "Résidus",
main = "Résidus vs ajustés")
# 3. Normalité des résidus
qqnorm(residuals(mod1))
qqline(residuals(mod1), col = "blue")
# 4. Autocorrélation des résidus (si données longitudinales)
acf(residuals(mod1), main = "ACF des résidus")
# 5. Normalité des pentes aléatoires (effets aléatoires sur time)
qqnorm(ranef(mod1)$ID[, "time"],
main = "QQ-plot des pentes aléatoires (time)")
qqline(ranef(mod1)$ID[, "time"], col = "blue")
#extrait P values
annot_pvals <- data.frame(
ADR1_time = "ADR1_J0",
value_log = max(F1$value_log, na.rm = TRUE) * 1.1,
label = "p[time] == 0.024 * '\n' * p[outcomes] == 0.066"
)
# ===========================
# Figure
# ===========================
Figure_1 <- ggplot(F1, aes(x = ADR1_time, y = value_log, fill = Outcomes)) +
geom_boxplot(position = position_dodge(0.7), width = 0.6) +
scale_y_log10() +
labs(
title = "Association between ADRB1 and outcomes",
x = "times of measurement",
y = "Percentage of T4 lymphocyte expressing ADRB1\n (log 10 scale)",
) +
geom_text(data = count_data, aes(x = ADR1_time, y =- min(F1$value_log, na.rm = TRUE) *0.1,
label = paste0("n=", n), group = Outcomes),
position = position_dodge(0.7), size = 4, vjust = 1)+
scale_x_discrete(labels = c("ADR1_J0" = "implantation", "ADR1_J3_J5" = "day 3 to 5", "ADR1_JS" = "explantation")) +
scale_fill_manual(values = c(
"ECMO Weaning" = "#4DAF4A",
"Bridge to transplant or LVAD" = "#377EB8",
"Death" = "#F8766D"
)) +
theme(legend.position = "none")+
annotate(
"text",
x = 1, # position sur l’axe x : 1 = "ADR1_J0"
y = max(F1$value_log, na.rm = TRUE) * 1.1, # position en haut du graphique
label = "p[time] == 0.024*';'~p[outcomes] == 0.066",
parse = TRUE,
size = 4,
hjust = 0
)
#####################
## ADRB2
#####################
F1_1 <- df %>%
select(ID,ADR2_J0, ADR2_J3_J5, ADR2_JS, Outcome)%>%
mutate(Outcomes = factor(case_when(Outcome=="Bridge to LVAD"|Outcome=="Bridge to Transplant"~ "Bridge to transplant or LVAD",T~Outcome), levels=c("Death","Bridge to transplant or LVAD","ECMO Weaning")))%>%
pivot_longer(
cols = c(ADR2_J0, ADR2_J3_J5, ADR2_JS),
names_to = "ADR2_time",
values_to = "value"
)%>%
mutate(time = as.numeric(case_when(ADR2_time=="ADR2_J0"~1,
ADR2_time=="ADR2_J3_J5"~2,
ADR2_time=="ADR2_JS"~3)),
value_log = log10(value))
count_data <- F1_1 %>%
filter(!is.na(value))%>%
group_by(ADR2_time, Outcomes) %>%
summarise(n = n(), .groups = "drop")%>%
filter(!is.na(Outcomes))
# ===========================
# Modélisation
# ===========================
# Modèle avec interaction
mod1 <- lmerTest::lmer(
value_log ~ time * Outcomes + (1 + time | ID),
data = F1_1
)
# Résumé du modèle
summary(mod1)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: value_log ~ time * Outcomes + (1 + time | ID)
Data: F1_1
REML criterion at convergence: 101.9
Scaled residuals:
Min 1Q Median 3Q Max
-2.25460 -0.42375 -0.02913 0.60839 1.65730
Random effects:
Groups Name Variance Std.Dev. Corr
ID (Intercept) 0.15994 0.3999
time 0.02624 0.1620 -0.77
Residual 0.08183 0.2861
Number of obs: 95, groups: ID, 46
Fixed effects:
Estimate Std. Error df t value
(Intercept) -0.02846 0.15731 49.17843 -0.181
time 0.08396 0.08554 43.13429 0.982
OutcomesBridge to transplant or LVAD 0.08591 0.29201 34.60684 0.294
OutcomesECMO Weaning 0.26694 0.20966 41.35055 1.273
time:OutcomesBridge to transplant or LVAD 0.12482 0.14021 28.85470 0.890
time:OutcomesECMO Weaning -0.04361 0.10803 35.41604 -0.404
Pr(>|t|)
(Intercept) 0.857
time 0.332
OutcomesBridge to transplant or LVAD 0.770
OutcomesECMO Weaning 0.210
time:OutcomesBridge to transplant or LVAD 0.381
time:OutcomesECMO Weaning 0.689
Correlation of Fixed Effects:
(Intr) time OBttoL OECMOW t:ttoL
time -0.868
OtcBttoLVAD -0.539 0.468
OtcmsECMOWn -0.750 0.651 0.404
t:OBttoLVAD 0.530 -0.610 -0.860 -0.397
tm:OtcECMOW 0.688 -0.792 -0.370 -0.858 0.483
# Modèle sans interaction (modèle réduit)
mod_no_interaction <- lmerTest::lmer(
value_log ~ time + Outcomes + (1 + time | ID),
data = F1_1
)
# Comparaison des modèles (test de l'interaction)
interact <- anova(mod_no_interaction, mod1)
# Interprétation : p = 0.76 → interaction non significative
# ===========================
# Tests des effets fixes
# ===========================
# ANOVA avec approximation de Kenward-Roger pour les p-values
a <- anova(mod_no_interaction, ddf = "Kenward-Roger")
print(a)
Type III Analysis of Variance Table with Kenward-Roger's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
time 0.25391 0.25391 1 33.452 3.0552 0.08965 .
Outcomes 0.45728 0.22864 2 36.789 2.7496 0.07711 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ===========================
# Vérification des hypothèses
# ===========================
# 1. Diagrammes de diagnostic généraux
plot(mod1) # résidus vs prédictions par groupe
# 2. Résidus vs valeurs ajustées (homoscédasticité)
plot(fitted(mod1), residuals(mod1),
xlab = "Valeurs ajustées", ylab = "Résidus",
main = "Résidus vs ajustés")
# 3. Normalité des résidus
qqnorm(residuals(mod1))
qqline(residuals(mod1), col = "blue")
# 4. Autocorrélation des résidus (si données longitudinales)
acf(residuals(mod1), main = "ACF des résidus")
# 5. Normalité des pentes aléatoires (effets aléatoires sur time)
qqnorm(ranef(mod1)$ID[, "time"],
main = "QQ-plot des pentes aléatoires (time)")
qqline(ranef(mod1)$ID[, "time"], col = "blue")
#extrait P values
annot_pvals <- data.frame(
ADR2_time = "ADR2_J0",
value_log = max(F1_1$value_log, na.rm = TRUE) * 1.1,
label = "p[time] == 0.07 * '\n' * p[outcomes] == 0.05"
)
Figure_1.1 <- ggplot(F1_1, aes(x = ADR2_time, y = value_log, fill = Outcomes)) +
geom_boxplot(position = position_dodge(0.7), width = 0.6) +
scale_y_log10() +
labs(
title = "Association between ADRB2 and outcomes",
x = "times of measurement",
y = "Percentage of T4 lymphocyte expressing ADRB2\n (log 10 scale)",
) +
geom_text(data = count_data, aes(x = ADR2_time, y = -min(F1_1$value_log, na.rm = TRUE) * 0.01,
label = paste0("n=", n), group = Outcomes),
position = position_dodge(0.7), size = 4, vjust = 1)+
scale_x_discrete(labels = c("ADR2_J0" = "implantation", "ADR2_J3_J5" = "day 3 to 5", "ADR2_JS" = "explantation"))+
scale_fill_manual(values = c(
"ECMO Weaning" = "#4DAF4A",
"Bridge to transplant or LVAD" = "#377EB8",
"Death" = "#F8766D"
)) +
annotate(
"text",
x = 1, # position sur l’axe x : 1 = "ADR1_J0"
y = max(F1_1$value_log, na.rm = TRUE) * 1.1, # position en haut du graphique
label = "p[time] == 0.07*';'~p[outcomes] == 0.05",
parse = TRUE,
size = 4,
hjust = 0
)
Fig_1 <- Figure_1 #+ theme(plot.title = element_blank())
Fig_1.1 <- Figure_1.1 #+ theme(plot.title = element_blank())
Figure1_tot <- Fig_1 + Fig_1.1 +
plot_annotation(title = "Association between percentage of T4 lymphocyte expressing ADRB and Outcome")
Figure1_tot
# Prepare data
F1m <- df %>%
select(ID, ADR1m_J0, ADR1m_J3_J5, ADR1m_JS, Outcome) %>%
mutate(Outcomes = factor(case_when(
Outcome == "Bridge to LVAD" | Outcome == "Bridge to Transplant" ~ "Bridge to transplant or LVAD",
TRUE ~ Outcome
), levels = c("Death", "Bridge to transplant or LVAD", "ECMO Weaning"))) %>%
pivot_longer(
cols = c(ADR1m_J0, ADR1m_J3_J5, ADR1m_JS),
names_to = "ADR1_time",
values_to = "value"
) %>%
mutate(
time = case_when(
ADR1_time == "ADR1m_J0" ~ 1,
ADR1_time == "ADR1m_J3_J5" ~ 2,
ADR1_time == "ADR1m_JS" ~ 3
)
) %>%
filter(!is.na(value), is.finite(value))
F1m <- F1m %>%
mutate(value_log = log10(value)) %>%
filter(!is.na(value_log), is.finite(value_log))
# Count data for annotation
count_data_m <- F1m %>%
filter(!is.na(value)) %>%
group_by(ADR1_time, Outcomes) %>%
summarise(n = n(), .groups = "drop") %>%
filter(!is.na(Outcomes))
# ===========================
# Mixed effects model
# ===========================
# Model with interaction
mod1_m <- lmerTest::lmer(
value_log ~ time * Outcomes + (1 + time | ID),
data = F1m
)
# Model without interaction
mod_no_interaction_m <- lmerTest::lmer(
value_log ~ time + Outcomes + (1 + time | ID),
data = F1m
)
# Interaction test
interact_m <- anova(mod_no_interaction_m, mod1_m)
# Interpretation : p = 0.64 → interaction not significant
# ANOVA for fixed effect (withouth interaction)
a_m <- anova(mod_no_interaction_m, ddf = "Kenward-Roger")
print(a_m)
Type III Analysis of Variance Table with Kenward-Roger's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
time 0.056058 0.056058 1 33.493 0.3582 0.5535
Outcomes 0.105353 0.052676 2 39.504 0.3365 0.7163
# ===========================
# Diagnostic plots
# ===========================
# 1. Residual diagnostic plots
plot(mod1_m)
# 2. Residuals vs Fitted
plot(fitted(mod1_m), residuals(mod1_m),
xlab = "Fitted values", ylab = "Residuals",
main = "Residuals vs Fitted")
# 3. Normality of residuals
qqnorm(residuals(mod1_m))
qqline(residuals(mod1_m), col = "blue")
# 4. Autocorrelation of residuals
acf(residuals(mod1_m), main = "ACF of residuals")
# 5. Normality of random intercepts (optional)
qqnorm(ranef(mod1_m)$ID[, "(Intercept)"],
main = "QQ plot of random intercepts (ID)")
qqline(ranef(mod1_m)$ID[, "(Intercept)"], col = "blue")
# ===========================
# Figure
# ===========================
Figure_1m <- ggplot(F1m, aes(
x = ADR1_time,
y = value,
fill = Outcomes,
group = interaction(ADR1_time, Outcomes)
)) +
geom_boxplot(position = position_dodge(0.7), width = 0.6) +
scale_y_log10() +
labs(
title = "Association between ADRB1 expression in monocytes and outcomes",
x = "Time of measurement",
y = "Percentage of monocytes expressing ADRB1\n(log10 scale)"
) +
geom_text(data = count_data_m, aes(
x = ADR1_time,
y = min(F1m$value[F1m$value > 0], na.rm = TRUE) * 0.5,
label = paste0("n=", n),
group = interaction(ADR1_time, Outcomes)
),
position = position_dodge(0.7),
size = 4,
vjust = 1) +
scale_x_discrete(labels = c(
"ADR1m_J0" = "implantation",
"ADR1m_J3_J5" = "day 3 to 5",
"ADR1m_JS" = "explantation"
)) +
scale_fill_manual(values = c(
"ECMO Weaning" = "#4DAF4A",
"Bridge to transplant or LVAD" = "#377EB8",
"Death" = "#F8766D"
)) +
theme(legend.position = "none") +
annotate(
"text",
x = 1,
y = max(F1m$value, na.rm = TRUE) * 1.1,
label = "p[time] == 0.955*';'~p[outcomes] == 0.649",
parse = TRUE,
size = 4,
hjust = 0
)
#####################
## ADRB2 (monocytes)
#####################
# ==========================
# Prepare data (non log)
# ==========================
F1m_2 <- df %>%
select(ID, ADR2m_J0, ADR2m_J3_J5, ADR2m_JS, Outcome) %>%
mutate(
Outcomes = factor(case_when(
Outcome %in% c("Bridge to LVAD", "Bridge to Transplant") ~ "Bridge to transplant or LVAD",
TRUE ~ Outcome
), levels = c("Death", "Bridge to transplant or LVAD", "ECMO Weaning"))
) %>%
pivot_longer(
cols = c(ADR2m_J0, ADR2m_J3_J5, ADR2m_JS),
names_to = "ADR2_time",
values_to = "value"
) %>%
mutate(
time = case_when(
ADR2_time == "ADR2m_J0" ~ 1,
ADR2_time == "ADR2m_J3_J5" ~ 2,
ADR2_time == "ADR2m_JS" ~ 3
)
) %>%
filter(!is.na(value), is.finite(value))
F1m_2 <- F1m_2 %>%
mutate(value_log = log10(value)) %>%
filter(!is.na(value_log), is.finite(value_log))
# ==========================
# Count available data
# ==========================
count_data_m2 <- F1m_2 %>%
group_by(ADR2_time, Outcomes) %>%
summarise(n = n(), .groups = "drop")
# ==========================
# Mixed effects model (simplified)
# ==========================
mod1_m2 <- lmerTest::lmer(
value_log ~ time * Outcomes + (1 + time | ID),
data = F1m_2
)
mod_no_interaction_m2 <- lmerTest::lmer(
value_log ~ time + Outcomes + (1 + time | ID),
data = F1m_2
)
# Test interaction
interact_m2 <- anova(mod_no_interaction_m2, mod1_m2)
# Interpretation : p = 0.76 → interaction not significant
# Fixed effect ANOVA (no interaction)
a_m2 <- anova(mod_no_interaction_m2, ddf = "Kenward-Roger")
print(a_m2)
Type III Analysis of Variance Table with Kenward-Roger's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
time 0.33421 0.33421 1 32.054 1.9455 0.1727
Outcomes 0.09096 0.04548 2 40.083 0.2647 0.7688
# Extract p-values
p_time <- a_m2["time", "Pr(>F)"]
p_outcomes <- a_m2["Outcomes", "Pr(>F)"]
annot_label <- paste0(
"p[time] == ", format(round(p_time, 3), nsmall = 3),
"*';'*~p[outcomes] == ", format(round(p_outcomes, 3), nsmall = 3)
)
# ==========================
# Diagnostics
# ==========================
plot(mod1_m2)
plot(fitted(mod1_m2), residuals(mod1_m2),
xlab = "Fitted values", ylab = "Residuals",
main = "Residuals vs Fitted")
qqnorm(residuals(mod1_m2))
qqline(residuals(mod1_m2), col = "blue")
acf(residuals(mod1_m2), main = "ACF of residuals")
qqnorm(ranef(mod1_m2)$ID[, "(Intercept)"], main = "QQ plot of random intercepts")
qqline(ranef(mod1_m2)$ID[, "(Intercept)"], col = "blue")
# ==========================
# Final Plot
# ==========================
Figure_1m_2 <- ggplot(F1m_2 %>% filter(value > 0), aes(
x = ADR2_time, y = value, fill = Outcomes,
group = interaction(ADR2_time, Outcomes)
)) +
geom_boxplot(position = position_dodge(0.7), width = 0.6) +
scale_y_log10() +
labs(
title = "Association between ADRB2 expression in monocytes and outcomes",
x = "Time of measurement",
y = "Percentage of monocytes expressing ADRB2\n(log10 scale)"
) +
geom_text(data = count_data_m2, aes(
x = ADR2_time,
y = min(F1m_2$value[F1m_2$value > 0], na.rm = TRUE) * 0.5,
label = paste0("n=", n),
group = interaction(ADR2_time, Outcomes)
),
position = position_dodge(0.7),
size = 4,
vjust = 1) +
scale_x_discrete(labels = c(
"ADR2m_J0" = "implantation",
"ADR2m_J3_J5" = "day 3 to 5",
"ADR2m_JS" = "explantation"
)) +
scale_fill_manual(values = c(
"ECMO Weaning" = "#4DAF4A",
"Bridge to transplant or LVAD" = "#377EB8",
"Death" = "#F8766D"
)) +
annotate(
"text",
x = 1,
y = max(F1m_2$value, na.rm = TRUE) * 1.2,
label = annot_label,
parse = TRUE,
size = 4,
hjust = 0
)
Fig_1m <- Figure_1m # + theme(plot.title = element_blank())
Fig_1m_2 <- Figure_1m_2 # + theme(plot.title = element_blank())
Figure1m_tot <- Fig_1m + Fig_1m_2 +
plot_annotation(title = "Association between percentage of monocytes expressing ADRB and Outcome")
Figure1m_tot
To explore the distribution of ADRB1 and ADRB2 expression over time in relation to clinical outcomes, boxplots were generated at three timepoints (implantation, day 3–5, and explantation). To evaluate the temporal dynamics of ADRB1 and ADRB2 expression across clinical outcomes, we used linear mixed-effects models with time, outcome group (ECMO weaning, bridge to LVAD/transplant, or death), and their interaction as fixed effects. Random intercepts and slopes for time were included to account for intra-patient variability. When appropriate, receptor expression values were log10-transformed. Interaction terms were tested using likelihood ratio tests between full and reduced models, and Type III ANOVA with Kenward-Roger approximation was used for fixed effects.
In patients who were successfully weaned from ECMO, both ADRB1 and ADRB2 expression in T4 lymphocytes showed an increase from implantation to explantation. In contrast, patients who died during ECMO exhibited lower or decreasing receptor expression over time. Although these trends were visually distinguishable, the interaction between time and outcome group was not statistically significant (ADRB1: p = 0.58; ADRB2: p = 0.76). Nevertheless, we observed a significant main effect of time for ADRB1 (p = 0.02) and a trend for ADRB2 (p = 0.07), alongside borderline differences across outcome groups (ADRB1: p = 0.07; ADRB2: p = 0.05).
In monocytes, ADRB1 and ADRB2 expression did not show any consistent time-dependent pattern across outcome groups. Boxplots revealed no evident increase or decrease over time, regardless of clinical trajecotry. These observations were supported by linear mixed-effects models, which showed no significant interaction between time and outcome group (ADRB1: pinteraction = 0.64; ADRB2: pinteraction = 0.76). Furthermore, neither time nor outcome group had significant main effects (ADRB1: p[time] = 0.95, p[outcomes] = 0.65; ADRB2: p[time] = 0.94, p[outcomes] = 0.57).
#ADRB1
F2_1_a <- df %>%
select(`jsev_epr_sev_t1_1,5l_itv`, ADR1_JS)%>%
mutate(
VTI_JS = case_when(
`jsev_epr_sev_t1_1,5l_itv` > 12 ~ ">12 cm",
`jsev_epr_sev_t1_1,5l_itv` <= 12 ~ "≤12 cm"
)
)%>%
na.omit()
# Mann-Whitney test for VTI and ADRB1
mann_whitney_VTI <- wilcox.test(ADR1_JS ~ VTI_JS, data = F2_1_a)
# Extract and format the p-value
p_value_VTI <- sprintf("%.3f", mann_whitney_VTI$p.value)
count_data <- F2_1_a %>%
group_by(VTI_JS) %>%
summarise(n = n())
# Boxplot for ADRB1 with p-value annotation
Figure2a <- ggplot(F2_1_a, aes(x = VTI_JS, y = log10(ADR1_JS), fill = VTI_JS)) +
geom_boxplot() +
geom_point(size = 2, alpha = 0.7, position = position_jitter(width = 0.1)) +
geom_text(data = count_data,
aes(x = VTI_JS, y = min(log10(F2_1_a$ADR1_JS), na.rm = TRUE) * 1.25,
label = paste0("n=", n)),
size = 4,
inherit.aes = FALSE) +
scale_fill_manual(values = c("≤12 cm" = "#0072B2", ">12 cm" = "#D55E00")) +
labs(
title = "A",
x = "VTI measured at 1.5 L/min",
y = "ADRB1 in lymphocyte measured at ECMO weaning \n (log10 scale)"
) +
theme_minimal() +
theme(legend.position = "none") +
annotate("text",
x = 1.5,
y = max(log10(F2_1_a$ADR1_JS), na.rm = TRUE),
label = paste("p:", p_value_VTI),
size = 5, color = "black", hjust = 0.5)
#ADRB2
F2_2_a <- df %>%
select(`jsev_epr_sev_t1_1,5l_itv`, ADR2_JS)%>%
mutate(
VTI_JS = case_when(
`jsev_epr_sev_t1_1,5l_itv` > 12 ~ ">12 cm",
`jsev_epr_sev_t1_1,5l_itv` <= 12 ~ "≤12 cm"
)
)%>%
na.omit()
# Mann-Whitney test for VTI and ADRB2
mann_whitney_VTI_2 <- wilcox.test(log10(ADR2_JS) ~ VTI_JS, data = F2_2_a)
# Extract and format the p-value
p_value_VTI_2 <- sprintf("%.3f", mann_whitney_VTI_2$p.value)
count_data <- F2_2_a %>%
group_by(VTI_JS) %>%
summarise(n = n(), .groups = "drop")
# Boxplot for ADRB2 with p-value annotation
Figure2b <- ggplot(F2_2_a, aes(x = VTI_JS, y = log10(ADR2_JS), fill = VTI_JS)) +
geom_boxplot() +
geom_point(size = 2, alpha = 0.7, position = position_jitter(width = 0.1)) +
scale_fill_manual(values = c("≤12 cm" = "#0072B2", ">12 cm" = "#D55E00")) +
labs(
title = "B",
x = "VTI measured at 1.5L /min",
y = "ADRB2 in lymphocyte measured at ECMO weaning \n (log10 scale)"
) +
theme_minimal() +
theme(legend.position = "none") +
annotate(
"text",
x = 1.5,
y = max(log10(F2_2_a$ADR2_JS), na.rm = TRUE),
label = paste("p:", p_value_VTI_2),
size = 5,
color = "black",
hjust = 0.5
) +
geom_text(
data = count_data,
aes(
x = VTI_JS,
y = min(log10(F2_2_a$ADR2_JS), na.rm = TRUE) * 1.25,
label = paste0("n=", n)
),
inherit.aes = FALSE,
size = 4
)
Figure2_tot <- Figure2a | Figure2b +
plot_annotation(
title = "Association between VTI at weaning and ADRB1/ADRB2 expression in lymphocyte",
tag_levels = c("A","B")
)
Figure2_tot
# ADRB1
F2_1_m <- df %>%
select(`jsev_epr_sev_t1_1,5l_itv`, ADR1m_JS) %>%
mutate(
VTI_JS = case_when(
`jsev_epr_sev_t1_1,5l_itv` > 12 ~ ">12 cm",
`jsev_epr_sev_t1_1,5l_itv` <= 12 ~ "≤12 cm"
),
VTI_JS = factor(VTI_JS, levels = c("≤12 cm", ">12 cm"))
) %>%
na.omit()
# Mann–Whitney test
mw_test_adr1m <- wilcox.test(log10(ADR1m_JS) ~ VTI_JS, data = F2_1_m)
# Extract p-value
p_value_adr1m <- sprintf("%.3f", mw_test_adr1m$p.value)
count_data_m <- F2_1_m %>%
group_by(VTI_JS) %>%
summarise(n = n())
# Plot
Figure2m_A <- ggplot(F2_1_m, aes(x = VTI_JS, y = log10(ADR1m_JS), fill = VTI_JS)) +
geom_boxplot() +
geom_point(position = position_jitter(width = 0.1), size = 2, alpha = 0.7) +
geom_text(data = count_data_m,
aes(x = VTI_JS,
y = -2.5,
label = paste0("n=", n)),
size = 4,
inherit.aes = FALSE) +
scale_fill_manual(values = c("≤12 cm" = "#0072B2", ">12 cm" = "#D55E00")) +
labs(
title = "A",
x = "VTI measured at 1.5L /min",
y = "ADRB1 in monocyte at ECMO weaning\n(log10 scale)"
) +
theme_minimal() +
theme(legend.position = "none") +
annotate(
"text",
x = 1.5,
y = max(log10(F2_1_m$ADR1m_JS), na.rm = TRUE) * 1.05,
label = paste("p:", p_value_adr1m),
size = 5
) +
coord_cartesian(clip = "off")
# ADRB2
F2_2_m <- df %>%
select(`jsev_epr_sev_t1_1,5l_itv`, ADR2m_JS) %>%
mutate(
VTI_JS = case_when(
`jsev_epr_sev_t1_1,5l_itv` > 12 ~ ">12 cm",
`jsev_epr_sev_t1_1,5l_itv` <= 12 ~ "≤12 cm"
),
VTI_JS = factor(VTI_JS, levels = c("≤12 cm", ">12 cm"))
) %>%
na.omit()
# Mann–Whitney test
mw_test_adr2m <- wilcox.test(log10(ADR2m_JS) ~ VTI_JS, data = F2_2_m)
# Extract p-value
p_value_adr2m <- sprintf("%.3f", mw_test_adr2m$p.value)
count_data_m2 <- F2_2_m %>%
group_by(VTI_JS) %>%
summarise(n = n())
# Plot
Figure2m_B <- ggplot(F2_2_m, aes(x = VTI_JS, y = log10(ADR2m_JS), fill = VTI_JS)) +
geom_boxplot() +
geom_point(position = position_jitter(width = 0.1), size = 2, alpha = 0.7) +
geom_text(data = count_data_m2,
aes(x = VTI_JS,
y = min(log10(F2_2_m$ADR2m_JS), na.rm = TRUE) * 1.25,
label = paste0("n=", n)),
size = 4,
inherit.aes = FALSE) +
scale_fill_manual(values = c("≤12 cm" = "#0072B2", ">12 cm" = "#D55E00")) +
labs(
title = "B",
x = "VTI measured at 1.5L /min",
y = "ADRB2 in monocyte at ECMO weaning\n(log10 scale)"
) +
theme_minimal() +
theme(legend.position = "none") +
annotate(
"text",
x = 1.5,
y = max(log10(F2_2_m$ADR2m_JS), na.rm = TRUE) * 1.05,
label = paste("p:", p_value_adr2m),
size = 5
)
# Combined plot
Figure2m_tot <- Figure2m_A | Figure2m_B +
plot_annotation(
title = "Association between VTI at weaning and ADRB1/ADRB2 expression in monocyte",
tag_levels = "A"
)
Figure2m_tot
To assess the association between ADRB1 and ADRB2 expression at the time of ECMO weaning and left ventricular outflow tract velocity-time integral (VTI) measured at 1.5 L/min ECMO flow, patients were stratified according to a VTI threshold of 12 cm. Mann–Whitney U tests were used to compare log-transformed ADRB expression between the two groups (≤12 cm vs >12 cm). Results were visualized using boxplots.
In lymphocytes, both ADRB1 and ADRB2 expression levels at the time of ECMO weaning were significantly lower in patients with a suboptimal VTI (≤12 cm), a threshold commonly used to assess readiness for explantation. This difference was statistically significant for both markers (ADRB1: p = 0.04; ADRB2: p = 0.04), suggesting a potential association between impaired hemodynamic recovery and reduced β-adrenergic receptor expression.
In contrast, no statistically significant differences in ADRB1 or ADRB2 expression in monocytes were observed between patients with VTI ≤12 cm and those with VTI >12 cm (ADRB1: p = 0.73; ADRB2: p = 0.79).
### Correlation between ADRB1 in Lymphocyte and total dobutamine dose
F3c_a <- df %>%
select(ADR1_JS, dose_tot_dobu) %>%
filter(!is.na(ADR1_JS), !is.na(dose_tot_dobu))
cor_test_adr1 <- cor.test(F3c_a$ADR1_JS, F3c_a$dose_tot_dobu, method = "spearman")
rho_adr1 <- sprintf("%.4f", cor_test_adr1$estimate)
pval_adr1 <- sprintf("p = %.3f", cor_test_adr1$p.value)
Fig_3c <- ggplot(F3c_a, aes(x = dose_tot_dobu, y = ADR1_JS)) +
geom_point(color = "blue", alpha = 0.6) +
geom_smooth(method = "lm", color = "red", se = TRUE) +
labs(title = "Correlation between dobutamine dose and ADRB1 (lymphocytes)",
x = "Dobutamine dose (mg)",
y = "ADRB1 expression on T helper lymphocytes") +
annotate("text",
x = min(F3c_a$dose_tot_dobu, na.rm = TRUE),
y = max(F3c_a$ADR1_JS, na.rm = TRUE),
label = paste("Spearman rho:", rho_adr1, "\n", pval_adr1),
hjust = 0, size = 5, color = "black") +
theme_minimal()
### Correlation between ADRB2 in Lymphocyte and total dobutamine dose
F3c_b <- df %>%
select(ADR2_JS, dose_tot_dobu) %>%
filter(!is.na(ADR2_JS), !is.na(dose_tot_dobu))
cor_test_adr2 <- cor.test(F3c_b$ADR2_JS, F3c_b$dose_tot_dobu, method = "spearman")
rho_adr2 <- sprintf("%.4f", cor_test_adr2$estimate)
pval_adr2 <- sprintf("p = %.3f", cor_test_adr2$p.value)
Fig_3d <- ggplot(F3c_b, aes(x = dose_tot_dobu, y = ADR2_JS)) +
geom_point(color = "blue", alpha = 0.6) +
geom_smooth(method = "lm", color = "red", se = TRUE) +
labs(title = "Correlation between dobutamine dose and ADRB2 (lymphocytes)",
x = "Dobutamine dose (mg)",
y = "ADRB2 expression on T helper lymphocytes") +
annotate("text",
x = min(F3c_b$dose_tot_dobu, na.rm = TRUE),
y = max(F3c_b$ADR2_JS, na.rm = TRUE),
label = paste("Spearman rho:", rho_adr2, "\n", pval_adr2),
hjust = 0, size = 5, color = "black") +
theme_minimal()
### Correlation between ADRB1 in Monocyte and total dobutamine dose
F3c_m1 <- df %>%
select(ADR1m_JS, dose_tot_dobu) %>%
filter(!is.na(ADR1m_JS), !is.na(dose_tot_dobu))
cor_test_adr1m <- cor.test(F3c_m1$ADR1m_JS, F3c_m1$dose_tot_dobu, method = "spearman")
rho_adr1m <- sprintf("%.4f", cor_test_adr1m$estimate)
pval_adr1m <- sprintf("p = %.3f", cor_test_adr1m$p.value)
Fig_3c_m <- ggplot(F3c_m1, aes(x = dose_tot_dobu, y = ADR1m_JS)) +
geom_point(color = "blue", alpha = 0.6) +
geom_smooth(method = "lm", color = "red", se = TRUE) +
labs(title = "Correlation between dobutamine dose and ADRB1 (monocytes)",
x = "Dobutamine dose (mg)",
y = "ADRB1 expression on monocytes at ECMO weaning") +
annotate("text",
x = min(F3c_m1$dose_tot_dobu, na.rm = TRUE),
y = max(F3c_m1$ADR1m_JS, na.rm = TRUE),
label = paste("Spearman rho:", rho_adr1m, "\n", pval_adr1m),
hjust = 0, size = 5, color = "black") +
theme_minimal()
### Correlation between ADRB2 in Monocyte and total dobutamine dose
F3c_m2 <- df %>%
select(ADR2m_JS, dose_tot_dobu) %>%
filter(!is.na(ADR2m_JS), !is.na(dose_tot_dobu))
cor_test_adr2m <- cor.test(F3c_m2$ADR2m_JS, F3c_m2$dose_tot_dobu, method = "spearman")
rho_adr2m <- sprintf("%.4f", cor_test_adr2m$estimate)
pval_adr2m <- sprintf("p = %.3f", cor_test_adr2m$p.value)
Fig_3d_m <- ggplot(F3c_m2, aes(x = dose_tot_dobu, y = ADR2m_JS)) +
geom_point(color = "blue", alpha = 0.6) +
geom_smooth(method = "lm", color = "red", se = TRUE) +
labs(title = "Correlation between dobutamine dose and ADRB2 (monocytes)",
x = "Dobutamine dose (mg)",
y = "ADRB2 expression on monocytes at ECMO weaning") +
annotate("text",
x = min(F3c_m2$dose_tot_dobu, na.rm = TRUE),
y = max(F3c_m2$ADR2m_JS, na.rm = TRUE),
label = paste("Spearman rho:", rho_adr2m, "\n", pval_adr2m),
hjust = 0, size = 5, color = "black") +
theme_minimal()
The association between cumulative dobutamine dose and the percentage of T-helper lymphocytes/monocyte expressing ADRB was assessed using Spearman’s rank correlation coefficient.
A weak negative correlation was observed between ADRB expression in lymphocytes and total dobutamine dose, which did not reach statistical significance. In contrast, both ADRB1 and ADRB2 expression in monocytes at the time of ECMO weaning were moderately correlated with cumulative dobutamine dose (Spearman ρ = 0.52, p = 0.021 for ADRB1; ρ = 0.54, p = 0.015 for ADRB2), indicating a significant positive association.
# ADRB1
F3_a <- df %>%
select(dose_tot_dobu, ADR1_JS) %>%
mutate(dobu_cut = case_when(
dose_tot_dobu > median(dose_tot_dobu, na.rm = TRUE) ~ "> median",
dose_tot_dobu <= median(dose_tot_dobu, na.rm = TRUE) ~ "≤ median"
)) %>%
filter(!is.na(ADR1_JS), !is.na(dose_tot_dobu))%>%
mutate(dobu_cut = factor(dobu_cut, levels = c("≤ median", "> median")))
# Test di Mann–Whitney
mw_adr1 <- wilcox.test(log10(ADR1_JS) ~ dobu_cut, data = F3_a)
p_val_adr1 <- sprintf("%.3f", mw_adr1$p.value)
count_data_a <- F3_a %>%
group_by(dobu_cut) %>%
summarise(n = n(), .groups = "drop")
# Figure ADR1
Figure3a <- ggplot(F3_a, aes(x = dobu_cut, y = log10(ADR1_JS), fill = dobu_cut)) +
geom_boxplot() +
geom_point(position = position_jitter(width = 0.1), alpha = 0.7, size = 2) +
geom_text(data = count_data_a,
aes(x = dobu_cut, y = min(log10(F3_a$ADR1_JS), na.rm = TRUE) * 1.1,
label = paste0("n=", n)),
size = 4,
inherit.aes = FALSE) +
scale_fill_manual(values = c("≤ median" = "#009E73", "> median" = "#CC79A7")) +
labs(
title = "A",
x = "Cumulative Dobutamine (cut at median)",
y = "ADRB1 at ECMO weaning (log10 scale)"
) +
theme_minimal() +
theme(legend.position = "none") +
annotate("text",
x = 1.5,
y = max(log10(F3_a$ADR1_JS), na.rm = TRUE) * 1.15,
label = paste("p:", p_val_adr1),
size = 5)
# ADRB2
F3_b <- df %>%
select(dose_tot_dobu, ADR2_JS) %>%
mutate(dobu_cut = case_when(
dose_tot_dobu > median(dose_tot_dobu, na.rm = TRUE) ~ "> median",
dose_tot_dobu <= median(dose_tot_dobu, na.rm = TRUE) ~ "≤ median"
)) %>%
filter(!is.na(ADR2_JS), !is.na(dose_tot_dobu))%>%
mutate(dobu_cut = factor(dobu_cut, levels = c("≤ median", "> median")))
# Test
mw_adr2 <- wilcox.test(log10(ADR2_JS) ~ dobu_cut, data = F3_b)
p_val_adr2 <- sprintf("%.3f", mw_adr2$p.value)
count_data_b <- F3_b %>%
group_by(dobu_cut) %>%
summarise(n = n(), .groups = "drop")
# Figure ADR2
Figure3b <- ggplot(F3_b, aes(x = dobu_cut, y = log10(ADR2_JS), fill = dobu_cut)) +
geom_boxplot() +
geom_point(position = position_jitter(width = 0.1), alpha = 0.7, size = 2) +
geom_text(data = count_data_b,
aes(x = dobu_cut, y = min(log10(F3_b$ADR2_JS), na.rm = TRUE) * 1.1,
label = paste0("n=", n)),
size = 4,
inherit.aes = FALSE) +
scale_fill_manual(values = c("≤ median" = "#009E73", "> median" = "#CC79A7")) +
labs(
title = "B",
x = "Cumulative Dobutamine (cut at median)",
y = "ADRB2 at ECMO weaning (log10 scale)"
) +
theme_minimal() +
theme(legend.position = "none") +
annotate("text",
x = 1.5,
y = max(log10(F3_b$ADR2_JS), na.rm = TRUE) * 1.15,
label = paste("p:", p_val_adr2),
size = 5)
# === Combine ===
Figure3_dobu_cum <- Figure3a + Figure3b +
plot_annotation(title = "Association between cumulative dobutamine and ADR expression in lymphocytes")
Figure3_dobu_cum
## ADRB1
# Prepare data
F3_m_a <- df %>%
select(dose_tot_dobu, ADR1m_JS) %>%
mutate(
dobu_cut = case_when(
dose_tot_dobu > median(dose_tot_dobu, na.rm = TRUE) ~ "> median",
TRUE ~ "≤ median"
)
) %>%
filter(!is.na(ADR1m_JS), !is.na(dose_tot_dobu)) %>%
mutate(dobu_cut = factor(dobu_cut, levels = c("≤ median", "> median")))
# Mann–Whitney test
mw_test_adr1m <- wilcox.test(log10(ADR1m_JS) ~ dobu_cut, data = F3_m_a)
p_value_adr1m <- sprintf("%.3f", mw_test_adr1m$p.value)
# Count per group
count_data_adr1m <- F3_m_a %>%
group_by(dobu_cut) %>%
summarise(n = n())
# Plot
Figure3_dobu_cum_a <- ggplot(F3_m_a, aes(x = dobu_cut, y = log10(ADR1m_JS), fill = dobu_cut)) +
geom_boxplot() +
geom_point(position = position_jitter(width = 0.1), alpha = 0.7, size = 2) +
geom_text(
data = F3_m_a %>% group_by(dobu_cut) %>% summarise(n = n()) %>%
mutate(y_n_label = -2.5), # y personalizzata direttamente qui
aes(x = dobu_cut, y = y_n_label, label = paste0("n=", n)),
inherit.aes = FALSE, size = 4
) +
scale_fill_manual(values = c("≤ median" = "#009E73", "> median" = "#CC79A7")) +
labs(title = "A", x = "Cumulative Dobutamine", y = "ADRB1 in monocytes at weaning (log10)") +
theme_minimal() +
theme(legend.position = "none") +
annotate("text",
x = 1.5,
y = max(log10(F3_m_a$ADR1m_JS), na.rm = TRUE),
label = paste("p:", p_value_adr1m),
size = 5)
##ADRB2
# Prepare data
F3_m_b <- df %>%
select(dose_tot_dobu, ADR2m_JS) %>%
mutate(
dobu_cut = case_when(
dose_tot_dobu > median(dose_tot_dobu, na.rm = TRUE) ~ "> median",
TRUE ~ "≤ median"
)
) %>%
filter(!is.na(ADR2m_JS), !is.na(dose_tot_dobu))%>%
mutate(dobu_cut = factor(dobu_cut, levels = c("≤ median", "> median")))
# Mann–Whitney test
mw_test_adr2m <- wilcox.test(log10(ADR2m_JS) ~ dobu_cut, data = F3_m_b)
p_value_adr2m <- sprintf("%.3f", mw_test_adr2m$p.value)
# Count per group
count_data_adr2m <- F3_m_b %>%
group_by(dobu_cut) %>%
summarise(n = n())
# Plot
Figure3_dobu_cum_b <- ggplot(F3_m_b, aes(x = dobu_cut, y = log10(ADR2m_JS), fill = dobu_cut)) +
geom_boxplot() +
geom_point(position = position_jitter(width = 0.1), alpha = 0.7, size = 2) +
geom_text(data = count_data_adr2m, aes(x = dobu_cut, y = min(log10(F3_m_b$ADR2m_JS), na.rm = TRUE) * 1.25,
label = paste0("n=", n)), inherit.aes = FALSE, size = 4) +
scale_fill_manual(values = c("≤ median" = "#009E73", "> median" = "#CC79A7")) +
labs(title = "B", x = "Cumulative Dobutamine", y = "ADRB2 in monocytes at weaning (log10)") +
theme_minimal() +
theme(legend.position = "none") +
annotate("text", x = 1.5, y = max(log10(F3_m_b$ADR2m_JS), na.rm = TRUE), label = paste("p:", p_value_adr2m), size = 5)
# Cumulative
Figure3m_dobu_cum <- Figure3_dobu_cum_a + Figure3_dobu_cum_b +
plot_annotation(title = "Association between cumulative dobutamine and ADR expression in monocytes")
Figure3m_dobu_cum
As evolution of some patients allowed an early weaning before D5, we considered that D3-D5 and D weaning are the same time point. To explore whether cumulative exposure to inotropic stimulation was associated with β-adrenergic receptor modulation over time, we investigated the relationship between the total cumulative dose of dobutamine (expressed in mg) and ADRB1/ADRB2 expression at ECMO weaning day. This timepoint was selected as it reflects the net effect of continuous pharmacological stimulation throughout the ECMO course, thereby allowing the assessment of potential receptor adaptation (e.g., up- or downregulation).
First, we performed a correlation analysis between cumulative dobutamine dose and ADRB1/ADRB2 expression using Spearman’s rank correlation coefficient, due to the non-normal distribution of the data. Second, patients were stratified according to the median of cumulative dobutamine dose. Expression levels of ADRB1 and ADRB2 (in both CD4+ lymphocytes and monocytes) at ECMO weaning were compared between patients above versus below the median using the Mann–Whitney U test. A two-sided p-value < 0.05 was considered statistically significant.
Among lymphocytes, there was no significant difference in ADRB expression between patients with higher versus lower cumulative dobutamine exposure. ADRB1 levels were comparable across groups, while a trend toward lower ADRB2 expression was observed in patients who received higher doses of dobutamine.
In contrast, monocyte expression patterns followed a more consistent trend. Both ADRB1 and ADRB2 expression levels tended to be higher in patients exposed to higher cumulative doses of dobutamine, although this association did not reach statistical significance.
### Correlation between ADRB1 in Lymphocyte and total noradrenalin dose
df4 <- df %>%
dplyr::select(dose_tot_NAd, ADR1_JS, Weight) %>%
filter(!is.na(ADR1_JS), !is.na(dose_tot_NAd)) %>%
mutate(NAD_per_kg = dose_tot_NAd / Weight)
# Correlation
correlation_adr1 <- sprintf("%.4f", cor(df4$ADR1_JS, df4$dose_tot_NAd, method = "spearman", use = "complete.obs"))
# Plot ADRB1
Figure_4_NADcum_A <- ggplot(df4, aes(x = dose_tot_NAd, y = ADR1_JS)) +
geom_point(color = "blue", alpha = 0.6) +
geom_smooth(method = "lm", color = "red", se = TRUE) +
labs(
title = "Correlation between noradrenaline dose and ADRB1 (lymphocytes)",
x = "Cumulative noradrenaline dose (mg)",
y = "ADRB1 expression at ECMO weaning (JS)"
) +
annotate(
"text",
x = min(df4$dose_tot_NAd, na.rm = TRUE),
y = max(df4$ADR1_JS, na.rm = TRUE) * 0.95,
label = paste("Spearman correlation:", correlation_adr1),
hjust = 0,
size = 5,
color = "black"
) +
theme_minimal()
print(Figure_4_NADcum_A)
### Correlation between ADRB2 in Lymphocyte and total noradrenalin dose
df4 <- df %>%
dplyr::select(dose_tot_NAd, ADR2_JS, Weight) %>%
filter(!is.na(ADR2_JS), !is.na(dose_tot_NAd)) %>%
mutate(NAD_per_kg = dose_tot_NAd / Weight)
correlation_adr2 <- sprintf("%.4f", cor(df4$ADR2_JS, df4$dose_tot_NAd, method = "spearman", use = "complete.obs"))
Figure_4_NADcum_B <- ggplot(df4, aes(x = dose_tot_NAd, y = ADR2_JS)) +
geom_point(color = "blue", alpha = 0.6) +
geom_smooth(method = "lm", color = "red", se = TRUE) +
labs(
title = "Correlation between noradrenaline dose and ADRB2 (lymphocytes)",
x = "Cumulative noradrenaline dose (mg)",
y = "ADRB2 expression at ECMO weaning (JS)"
) +
annotate(
"text",
x = min(df4$dose_tot_NAd, na.rm = TRUE),
y = max(df4$ADR2_JS, na.rm = TRUE) * 0.95,
label = paste("Spearman correlation:", correlation_adr2),
hjust = 0,
size = 5,
color = "black"
) +
theme_minimal()
print(Figure_4_NADcum_B)
### Correlation between ADRB1 in Monocyte and total noradrenalin dose
df4 <- df %>%
select(ADR1m_JS, dose_tot_NAd, Weight) %>%
filter(!is.na(ADR1m_JS), !is.na(dose_tot_NAd)) %>%
mutate(NAD_per_kg = dose_tot_NAd / Weight)
# Spearman correlation
cor_adr1m <- sprintf("%.4f", cor(df4$ADR1m_JS, df4$dose_tot_NAd, method = "spearman", use = "complete.obs"))
# Plot
Figure4_NADmono_A <- ggplot(df4, aes(x = dose_tot_NAd, y = ADR1m_JS)) +
geom_point(color = "blue", alpha = 0.6) +
geom_smooth(method = "lm", color = "red", se = TRUE) +
labs(
title = "Correlation between Noradrenaline dose and ADRB1 (monocytes)",
x = "Cumulative noradrenaline dose (mg)",
y = "Monocyte ADRB1 expression at ECMO weaning"
) +
annotate("text",
x = min(df4$dose_tot_NAd, na.rm = TRUE),
y = max(df4$ADR1m_JS, na.rm = TRUE) * 0.95,
label = paste("Spearman correlation:", cor_adr1m),
hjust = 0, size = 5, color = "black") +
theme_minimal()
print(Figure4_NADmono_A)
### Correlation between ADRB2 in Monocyte and total noradrenalin dose
df4 <- df %>%
select(ADR2m_JS, dose_tot_NAd, Weight) %>%
filter(!is.na(ADR2m_JS), !is.na(dose_tot_NAd)) %>%
mutate(NAD_per_kg = dose_tot_NAd / Weight)
# Spearman correlation
cor_adr2m <- sprintf("%.4f", cor(df4$ADR2m_JS, df4$dose_tot_NAd, method = "spearman", use = "complete.obs"))
# Plot
Figure4_NADmono_B <- ggplot(df4, aes(x = dose_tot_NAd, y = ADR2m_JS)) +
geom_point(color = "blue", alpha = 0.6) +
geom_smooth(method = "lm", color = "red", se = TRUE) +
labs(
title = "Correlation between Noradrenaline dose and ADRB2 (monocytes)",
x = "Cumulative noradrenaline dose (mg)",
y = "Monocyte ADRB2 expression at ECMO weaning"
) +
annotate("text",
x = min(df4$dose_tot_NAd, na.rm = TRUE),
y = max(df4$ADR2m_JS, na.rm = TRUE) * 0.95,
label = paste("Spearman correlation:", cor_adr2m),
hjust = 0, size = 5, color = "black") +
theme_minimal()
print(Figure4_NADmono_B)
No significant correlations were observed between cumulative noradrenaline dose and ADRB1 or ADRB2 expression in either T helper lymphocytes (r = –0.09 and –0.29, respectively) or monocytes (r = –0.18 and 0.02, respectively).
### ADRB1 in T lymphocytes ~ cumulative NAd (median)
F5a <- df %>%
mutate(
ADR1_JS_clean = case_when(
!is.na(ADR1_JS) ~ ADR1_JS,
is.na(ADR1_JS) & !is.na(ADR1_J3_J5) ~ ADR1_J3_J5,
TRUE ~ NA_real_
),
NAd_group = case_when(
dose_tot_NAd > median(dose_tot_NAd, na.rm = TRUE) ~ "> median",
dose_tot_NAd <= median(dose_tot_NAd, na.rm = TRUE) ~ "≤ median",
TRUE ~ NA_character_
)
) %>%
filter(!is.na(ADR1_JS_clean), !is.na(NAd_group)) %>%
mutate(NAd_group = factor(NAd_group, levels = c("≤ median", "> median")))
# p-value (Mann–Whitney)
wilcox_test_adr1 <- wilcox.test(ADR1_JS_clean ~ NAd_group, data = F5a)
p_value_adr1 <- sprintf("%.3f", wilcox_test_adr1$p.value)
count_data_adr1 <- F5a %>%
group_by(NAd_group) %>%
summarise(n = n())
# Boxplot
Figure_5a <- ggplot(F5a, aes(x = NAd_group, y = ADR1_JS_clean, fill = NAd_group)) +
geom_boxplot() +
geom_point(position = position_jitter(width = 0.1), size = 2, alpha = 0.7) +
geom_text(
data = count_data_adr1,
aes(x = NAd_group, y = min(F5a$ADR1_JS_clean, na.rm = TRUE) * 1.2,
label = paste0("n=", n)),
size = 4,
inherit.aes = FALSE
) +
scale_fill_manual(values = c("≤ median" = "#66C2AC", "> median" = "#D73027")) +
labs(
title = "A",
x = "Cumulative noradrenaline dose",
y = "ADRB1 expression at ECMO weaning"
) +
annotate("text",
x = 1.5,
y = max(F5a$ADR1_JS_clean, na.rm = TRUE) * 1.05,
label = paste("p:", p_value_adr1),
size = 5, color = "black") +
theme_minimal() +
theme(legend.position = "none")
### ADRB2 in T lymphocytes ~ cumulative NAd (median)
F5b <- df %>%
mutate(
ADR2_JS_clean = case_when(
!is.na(ADR2_JS) ~ ADR2_JS,
is.na(ADR2_JS) & !is.na(ADR2_J3_J5) ~ ADR2_J3_J5,
TRUE ~ NA_real_
),
NAd_group = case_when(
dose_tot_NAd > median(dose_tot_NAd, na.rm = TRUE) ~ "> median",
dose_tot_NAd <= median(dose_tot_NAd, na.rm = TRUE) ~ "≤ median",
TRUE ~ NA_character_
)
) %>%
filter(!is.na(ADR2_JS_clean), !is.na(NAd_group)) %>%
mutate(NAd_group = factor(NAd_group, levels = c("≤ median", "> median")))
# Test di Mann–Whitney
wilcox_test_adr2 <- wilcox.test(ADR2_JS_clean ~ NAd_group, data = F5b)
p_value_adr2 <- sprintf("%.3f", wilcox_test_adr2$p.value)
# Conta per gruppo
count_data_adr2 <- F5b %>%
group_by(NAd_group) %>%
summarise(n = n())
# Boxplot
Figure_5b <- ggplot(F5b, aes(x = NAd_group, y = ADR2_JS_clean, fill = NAd_group)) +
geom_boxplot() +
geom_point(position = position_jitter(width = 0.1), size = 2, alpha = 0.7) +
geom_text(
data = count_data_adr2,
aes(x = NAd_group, y = min(F5b$ADR2_JS_clean, na.rm = TRUE) * 1.2,
label = paste0("n=", n)),
size = 4,
inherit.aes = FALSE
) +
scale_fill_manual(values = c("≤ median" = "#66C2AC", "> median" = "#D73027")) +
labs(
title = "B",
x = "Cumulative noradrenaline dose",
y = "ADRB2 expression at ECMO weaning"
) +
annotate("text",
x = 1.5,
y = max(F5b$ADR2_JS_clean, na.rm = TRUE) * 1.05,
label = paste("p:", p_value_adr2),
size = 5, color = "black") +
theme_minimal() +
theme(legend.position = "none")
Figure_5a <- Figure_5a + theme(plot.title = element_blank())
Figure_5b <- Figure_5b + theme(plot.title = element_blank())
# Combinazione dei due plot
Figure4_tot <- Figure_5a + Figure_5b +
plot_annotation(
title = "Association between cumulative noradrenaline dose and β-adrenergic receptor expression on T helper lymphocytes"
)
Figure4_tot
cut_median_nad <- median(df$dose_tot_NAd, na.rm = TRUE)
#ADRB1m
F4_m_a <- df %>%
select(dose_tot_NAd, ADR1m_JS, ADR1m_J3_J5) %>%
mutate(
ADR1m_final = case_when(
!is.na(ADR1m_JS) ~ ADR1m_JS,
is.na(ADR1m_JS) & !is.na(ADR1m_J3_J5) ~ ADR1m_J3_J5,
TRUE ~ NA_real_
),
NAd_cut = case_when(
dose_tot_NAd <= cut_median_nad ~ "≤ median",
dose_tot_NAd > cut_median_nad ~ "> median",
TRUE ~ NA_character_
)
) %>%
filter(!is.na(ADR1m_final), !is.na(NAd_cut)) %>%
mutate(NAd_cut = factor(NAd_cut, levels = c("≤ median", "> median")))
# P-value
p_val_adr1m <- wilcox.test(ADR1m_final ~ NAd_cut, data = F4_m_a)$p.value
p_text_adr1m <- paste("p:", format(round(p_val_adr1m, 3), nsmall = 3))
# Count
count_data_m1 <- F4_m_a %>%
group_by(NAd_cut) %>%
summarise(n = n())
# Plot ADRB1m
Figure4m_A <- ggplot(F4_m_a, aes(x = NAd_cut, y = log10(ADR1m_final), fill = NAd_cut)) +
geom_boxplot() +
geom_point(position = position_jitter(width = 0.1), size = 2, alpha = 0.7) +
geom_text(data = count_data_m1,
aes(x = NAd_cut, y = min(log10(F4_m_a$ADR1m_final), na.rm = TRUE) * 1.2,
label = paste0("n=", n)),
size = 4, inherit.aes = FALSE) +
scale_fill_manual(values = c("≤ median" = "#66C2A5", "> median" = "#D73027")) +
labs(
title = "A",
x = "Cumulative norepinephrine dose",
y = "Monocyte ADRB1 expression at ECMO weaning\n(log10 scale)"
) +
theme_minimal() +
theme(legend.position = "none") +
annotate("text", x = 1.5,
y = max(log10(F4_m_a$ADR1m_final), na.rm = TRUE) * 1.15,
label = p_text_adr1m,
size = 5)
# ADRB2m
F4_m_b <- df %>%
select(dose_tot_NAd, ADR2m_JS, ADR2m_J3_J5) %>%
mutate(
ADR2m_final = case_when(
!is.na(ADR2m_JS) ~ ADR2m_JS,
is.na(ADR2m_JS) & !is.na(ADR2m_J3_J5) ~ ADR2m_J3_J5,
TRUE ~ NA_real_
),
NAd_cut = case_when(
dose_tot_NAd <= cut_median_nad ~ "≤ median",
dose_tot_NAd > cut_median_nad ~ "> median",
TRUE ~ NA_character_
)
) %>%
filter(!is.na(ADR2m_final), !is.na(NAd_cut)) %>%
mutate(NAd_cut = factor(NAd_cut, levels = c("≤ median", "> median")))
# P-value
p_val_adr2m <- wilcox.test(ADR2m_final ~ NAd_cut, data = F4_m_b)$p.value
p_text_adr2m <- paste("p:", format(round(p_val_adr2m, 3), nsmall = 3))
# Count
count_data_m2 <- F4_m_b %>%
group_by(NAd_cut) %>%
summarise(n = n())
# Plot ADRB2m
Figure4m_B <- ggplot(F4_m_b, aes(x = NAd_cut, y = log10(ADR2m_final), fill = NAd_cut)) +
geom_boxplot() +
geom_point(position = position_jitter(width = 0.1), size = 2, alpha = 0.7) +
geom_text(data = count_data_m2,
aes(x = NAd_cut, y = min(log10(F4_m_b$ADR2m_final), na.rm = TRUE) * 1.2,
label = paste0("n=", n)),
size = 4, inherit.aes = FALSE) +
scale_fill_manual(values = c("≤ median" = "#66C2AC", "> median" = "#D73027")) +
labs(
title = "B",
x = "Cumulative norepinephrine dose",
y = "Monocyte ADRB2 expression at ECMO weaning\n(log10 scale)"
) +
theme_minimal() +
theme(legend.position = "none") +
annotate("text", x = 1.5,
y = max(log10(F4_m_b$ADR2m_final), na.rm = TRUE) * 1.15,
label = p_text_adr2m,
size = 5)
# Combine plots
Figure4_mono <- Figure4m_A + Figure4m_B +
plot_annotation(title = "Association between cumulative norepinephrine dose and ADRB1/2 expression in monocytes at ECMO weaning")
Figure4_mono
Among lymphocytes, both ADRB1 and ADRB2 expression at ECMO weaning tended to be lower in patients receiving a cumulative noradrenaline dose above the median. However, these differences did not reach statistical significance.
Among monocytes, ADRB1 expression at ECMO weaning was significantly lower in patients receiving a cumulative norepinephrine dose above the median compared to those below (p = 0.039). A similar trend was observed for ADRB2 expression, although the difference did not reach statistical significance (p = 0.077).
While Spearman correlation analyses did not reveal a significant association between cumulative norepinephrine dose and ADRB1/ADRB2 expression in monocytes, patients with doses above the median exhibited significantly lower ADRB1 levels at ECMO weaning (p = 0.039). This suggests that potential receptor downregulation may occur only above a certain threshold of catecholaminergic stimulation.
df <- df %>%
mutate(VIS_cut = case_when(
VIS_D0 > 15 ~ ">15",
VIS_D0 <= 15 ~ "≤15"
)) %>%
mutate(VIS_cut = factor(VIS_cut, levels = c("≤15", ">15")))
# ADRB1 in lymphocytes
F_VIS_adrb1 <- df %>%
select(ADR1_JS, VIS_cut) %>%
filter(!is.na(ADR1_JS), !is.na(VIS_cut))
count_data_adrb1 <- F_VIS_adrb1 %>%
group_by(VIS_cut) %>%
summarise(n = n(), .groups = "drop")
p_value_adrb1 <- wilcox.test(log10(ADR1_JS) ~ VIS_cut, data = F_VIS_adrb1)$p.value
p_label_adrb1 <- paste("p =", format(p_value_adrb1, digits = 2, nsmall = 3))
plot_ADRB1 <- ggplot(F_VIS_adrb1, aes(x = VIS_cut, y = log10(ADR1_JS), fill = VIS_cut)) +
geom_boxplot() +
geom_point(position = position_jitter(width = 0.1), size = 2, alpha = 0.7) +
geom_text(data = count_data_adrb1, aes(x = VIS_cut, y = -2.5, label = paste0("n=", n)),
inherit.aes = FALSE, size = 4) +
annotate("text", x = 1.5, y = max(log10(F_VIS_adrb1$ADR1_JS), na.rm = TRUE), label = p_label_adrb1, size = 5) +
scale_fill_manual(values = c("≤15" = "#0072B2", ">15" = "#D55E00")) +
labs(title = "ADRB1 expression in lymphocytes by VIS score", x = "VIS score", y = "ADRB1 in lymphocytes at day 0 (log10 scale)") +
theme_minimal() +
theme(legend.position = "none")
# ADRB2 in lymphocytes
F_VIS_adrb2 <- df %>%
select(ADR2_JS, VIS_cut) %>%
filter(!is.na(ADR2_JS), !is.na(VIS_cut))
count_data_adrb2 <- F_VIS_adrb2 %>%
group_by(VIS_cut) %>%
summarise(n = n(), .groups = "drop")
p_value_adrb2 <- wilcox.test(log10(ADR2_JS) ~ VIS_cut, data = F_VIS_adrb2)$p.value
p_label_adrb2 <- paste("p =", format(p_value_adrb2, digits = 2, nsmall = 3))
plot_ADRB2 <- ggplot(F_VIS_adrb2, aes(x = VIS_cut, y = log10(ADR2_JS), fill = VIS_cut)) +
geom_boxplot() +
geom_point(position = position_jitter(width = 0.1), size = 2, alpha = 0.7) +
geom_text(data = count_data_adrb2, aes(x = VIS_cut, y = -2.5, label = paste0("n=", n)),
inherit.aes = FALSE, size = 4) +
annotate("text", x = 1.5, y = max(log10(F_VIS_adrb2$ADR2_JS), na.rm = TRUE), label = p_label_adrb2, size = 5) +
scale_fill_manual(values = c("≤15" = "#0072B2", ">15" = "#D55E00")) +
labs(title = "ADRB2 expression in lymphocytes by VIS score", x = "VIS score", y = "ADRB2 in lymphocytes at day 0 (log10 scale)") +
theme_minimal() +
theme(legend.position = "none")
# Combined plot per linfociti
plot_ADRB_lympho_combined <- plot_ADRB1 + plot_ADRB2 +
plot_layout(ncol = 2) +
plot_annotation(title = "ADRB1 and ADRB2 expression in lymphocytes stratified by VIS score")
print(plot_ADRB_lympho_combined)
# VIS_cut
df <- df %>%
mutate(VIS_cut = case_when(
VIS_D0 > 15 ~ ">15",
VIS_D0 <= 15 ~ "≤15"
)) %>%
mutate(VIS_cut = factor(VIS_cut, levels = c("≤15", ">15")))
### ADRB1
F_VIS_adrb1_m <- df %>%
select(ADR1m_JS, VIS_cut) %>%
filter(!is.na(ADR1m_JS), !is.na(VIS_cut))
count_data_adr1m <- F_VIS_adrb1_m %>%
group_by(VIS_cut) %>%
summarise(n = n(), .groups = "drop")
p_value_adr1m <- wilcox.test(log10(ADR1m_JS) ~ VIS_cut, data = F_VIS_adrb1_m)$p.value
p_label_adr1m <- paste("p =", format(p_value_adr1m, digits = 2, nsmall = 3))
plot_ADRB1m <- ggplot(F_VIS_adrb1_m, aes(x = VIS_cut, y = log10(ADR1m_JS), fill = VIS_cut)) +
geom_boxplot() +
geom_point(position = position_jitter(width = 0.1), size = 2, alpha = 0.7) +
geom_text(data = count_data_adr1m, aes(x = VIS_cut, y = -2.5, label = paste0("n=", n)),
inherit.aes = FALSE, size = 4) +
annotate("text", x = 1.5, y = max(log10(F_VIS_adrb1_m$ADR1m_JS), na.rm = TRUE), label = p_label_adr1m, size = 5) +
scale_fill_manual(values = c("≤15" = "#0072B2", ">15" = "#D55E00")) +
labs(title = "ADRB1 expression in monocytes by VIS", x = "VIS score", y = "ADRB1 in monocytes at day 0 (log10 scale)") +
theme_minimal() +
theme(legend.position = "none")
### ADRB2
F_VIS_adrb2_m <- df %>%
select(ADR2m_JS, VIS_cut) %>%
filter(!is.na(ADR2m_JS), !is.na(VIS_cut))
count_data_adr2m <- F_VIS_adrb2_m %>%
group_by(VIS_cut) %>%
summarise(n = n(), .groups = "drop")
p_value_adr2m <- wilcox.test(log10(ADR2m_JS) ~ VIS_cut, data = F_VIS_adrb2_m)$p.value
p_label_adr2m <- paste("p =", format(p_value_adr2m, digits = 2, nsmall = 3))
plot_ADRB2m <- ggplot(F_VIS_adrb2_m, aes(x = VIS_cut, y = log10(ADR2m_JS), fill = VIS_cut)) +
geom_boxplot() +
geom_point(position = position_jitter(width = 0.1), size = 2, alpha = 0.7) +
geom_text(data = count_data_adr2m, aes(x = VIS_cut, y = -2.5, label = paste0("n=", n)),
inherit.aes = FALSE, size = 4) +
annotate("text", x = 1.5, y = max(log10(F_VIS_adrb2_m$ADR2m_JS), na.rm = TRUE), label = p_label_adr2m, size = 5) +
scale_fill_manual(values = c("≤15" = "#0072B2", ">15" = "#D55E00")) +
labs(title = "ADRB2 expression in monocytes by VIS", x = "VIS score", y = "ADRB2 in monocytes at day 0 (log10 scale)") +
theme_minimal() +
theme(legend.position = "none")
# Combined
plot_ADRB_mono_combined <- plot_ADRB1m + plot_ADRB2m +
plot_layout(ncol = 2) +
plot_annotation(title = "ADRB1 and ADRB2 expression in monocytes stratified by VIS score")
# Combined plot per linfociti
plot_ADRB_lympho_combined <- plot_ADRB1 + plot_ADRB2 +
plot_layout(ncol = 2) +
plot_annotation(title = "ADRB1 and ADRB2 expression in lymphocytes stratified by VIS score")
print(plot_ADRB_mono_combined)
#ADRB1
F5_lym <- df %>%
dplyr::select(Levosimendan, ADR1_J3_J5) %>%
filter(!is.na(Levosimendan), !is.na(ADR1_J3_J5)) %>%
mutate(log_ADR1 = log10(ADR1_J3_J5))
# p-value
p_adr1 <- sprintf("%.4f", wilcox.test(log_ADR1 ~ Levosimendan, data = F5_lym, exact = FALSE)$p.value)
# count
count_data_adr1 <- F5_lym %>%
group_by(Levosimendan) %>%
summarise(n = n(), .groups = "drop")
# Boxplot
Figure_5A <- ggplot(F5_lym, aes(x = Levosimendan, y = log_ADR1, fill = Levosimendan)) +
geom_boxplot() +
geom_point(position = position_jitter(width = 0.1), size = 2, alpha = 0.7) +
geom_text(data = count_data_adr1,
aes(x = Levosimendan,
y = min(F5_lym$log_ADR1, na.rm = TRUE) * 1.25,
label = paste0("n=", n)),
inherit.aes = FALSE,
size = 4) +
scale_fill_manual(values = c("NO" = "#56B4E9", "YES" = "#F0E442")) +
labs(title = "ADRB1 expression on T helper lymphocytes",
x = "Levosimendan",
y = "ADRB1 at D3–D5 of ECMO (log10 scale)") +
annotate("text", x = 1.5, y = max(F5_lym$log_ADR1, na.rm = TRUE),
label = paste("Mann–Whitney p =", p_adr1),
size = 5, hjust = 0.5, color = "black") +
theme_minimal() +
theme(legend.position = "none")
#ADRB2
F5_lym_b <- df %>%
dplyr::select(Levosimendan, ADR2_J3_J5) %>%
filter(!is.na(Levosimendan), !is.na(ADR2_J3_J5)) %>%
mutate(log_ADR2 = log10(ADR2_J3_J5))
# p-value
p_adr2 <- sprintf("%.4f", wilcox.test(log_ADR2 ~ Levosimendan, data = F5_lym_b, exact = FALSE)$p.value)
# count
count_data_adr2 <- F5_lym_b %>%
group_by(Levosimendan) %>%
summarise(n = n(), .groups = "drop")
# Boxplot
Figure_5B <- ggplot(F5_lym_b, aes(x = Levosimendan, y = log_ADR2, fill = Levosimendan)) +
geom_boxplot() +
geom_point(position = position_jitter(width = 0.1), size = 2, alpha = 0.7) +
geom_text(data = count_data_adr2,
aes(x = Levosimendan,
y = min(F5_lym_b$log_ADR2, na.rm = TRUE) * 1.25,
label = paste0("n=", n)),
inherit.aes = FALSE,
size = 4) +
scale_fill_manual(values = c("NO" = "#56B4E9", "YES" = "#F0E442")) +
labs(title = "ADRB2 expression on T helper lymphocytes",
x = "Levosimendan",
y = "ADRB2 at D3–D5 of ECMO (log10 scale)") +
annotate("text", x = 1.5, y = max(F5_lym_b$log_ADR2, na.rm = TRUE),
label = paste("Mann–Whitney p =", p_adr2),
size = 5, hjust = 0.5, color = "black") +
theme_minimal() +
theme(legend.position = "none")
Figure_5A <- Figure_5A + theme(plot.title = element_blank())
Figure_5B <- Figure_5B + theme(plot.title = element_blank())
Figure5_tot_log <- Figure_5A + Figure_5B +
plot_annotation(title = "Association between Levosimendan and ADRB expression (log scale) at D3–D5")
Figure5_tot_log
# ADRB1
F5_mon_a <- df %>%
dplyr::select(Levosimendan, ADR1m_J3_J5) %>%
filter(!is.na(Levosimendan), !is.na(ADR1m_J3_J5))
# P-value
p_adr1m <- sprintf("%.4f", wilcox.test(ADR1m_J3_J5 ~ Levosimendan, data = F5_mon_a, exact = FALSE)$p.value)
# Count per group
count_data_adr1m <- F5_mon_a %>%
group_by(Levosimendan) %>%
summarise(n = n())
# Boxplot
Figure_5Am <- ggplot(F5_mon_a, aes(x = Levosimendan, y = log10(ADR1m_J3_J5), fill = Levosimendan)) +
geom_boxplot() +
geom_point(size = 2, alpha = 0.7, position = position_jitter(width = 0.1)) +
geom_text(data = count_data_adr1m,
aes(x = Levosimendan, y = min(log10(F5_mon_a$ADR1m_J3_J5), na.rm = TRUE) * 1.2,
label = paste0("n=", n)),
inherit.aes = FALSE,
size = 4) +
scale_fill_manual(values = c("NO" = "#56B4E9", "YES" = "#F0E442")) +
labs(title = "A",
x = "Levosimendan",
y = "Monocyte ADRB1 expression at D3–D5 (log10 scale)") +
annotate("text", x = 1.5, y = max(log10(F5_mon_a$ADR1m_J3_J5), na.rm = TRUE),
label = paste("Mann–Whitney p =", p_adr1m),
size = 5, color = "black", hjust = 0.5) +
theme_minimal() +
theme(legend.position = "none")
# ADRB2
F5_mon_b <- df %>%
dplyr::select(Levosimendan, ADR2m_J3_J5) %>%
filter(!is.na(Levosimendan), !is.na(ADR2m_J3_J5))
# P-value
p_adr2m <- sprintf("%.4f", wilcox.test(ADR2m_J3_J5 ~ Levosimendan, data = F5_mon_b, exact = FALSE)$p.value)
count_data_adr2m <- F5_mon_b %>%
group_by(Levosimendan) %>%
summarise(n = n())
# Boxplot
Figure_5Bm <- ggplot(F5_mon_b, aes(x = Levosimendan, y = log10(ADR2m_J3_J5), fill = Levosimendan)) +
geom_boxplot() +
geom_point(size = 2, alpha = 0.7, position = position_jitter(width = 0.1)) +
geom_text(data = count_data_adr2m,
aes(x = Levosimendan, y = min(log10(F5_mon_b$ADR2m_J3_J5), na.rm = TRUE) * 1.2,
label = paste0("n=", n)),
inherit.aes = FALSE,
size = 4) +
scale_fill_manual(values = c("NO" = "#56B4E9", "YES" = "#F0E442")) +
labs(title = "B",
x = "Levosimendan",
y = "Monocyte ADRB2 expression at D3–D5 (log10 scale)") +
annotate("text", x = 1.5, y = max(log10(F5_mon_b$ADR2m_J3_J5), na.rm = TRUE),
label = paste("Mann–Whitney p =", p_adr2m),
size = 5, color = "black", hjust = 0.5) +
theme_minimal() +
theme(legend.position = "none")
#Combined
Figure_5Am <- Figure_5Am + theme(plot.title = element_blank())
Figure_5Bm <- Figure_5Bm + theme(plot.title = element_blank())
Figure5_mon_tot <- Figure_5Am + Figure_5Bm +
plot_annotation(title = "Association between Levosimendan and ADRB expression on monocytes at D3–D5 of ECMO",
tag_levels = "A")
Figure5_mon_tot
ADRB2 expression at D3-D5 was significantly lower in patients who did not received Levosimendan compared to those who did. N difference where seen in monocytes.
# ADRB1
km_fit <- survfit(Surv(diff_days_J28, Outcome_censored_28) ~ ADR1_D0_median, data = df)
Figure_6 <- ggsurvplot(km_fit, data = df,
pval = TRUE,
pval.coord = c(10, 0.1),
conf.int = FALSE,
risk.table = TRUE,
risk.table.height = 0.2,
risk.table.y.text = TRUE,
risk.table.title = "Number at Risk (number censored)",
xlim = c(0, 28),
break.time.by = 7,
legend.title = " ",
palette = c("#E63946", "#457B9D"),
ggtheme = theme_minimal(),
legend.labs = c("High", "Low")
)
Figure_6 <- Figure_6 + ggtitle("28-day Survival Curves for ADRB1")
# ADRB2
km_fit <- survfit(Surv(diff_days_J28, Outcome_censored_28) ~ ADR2_D0_median, data = df)
Figure_6.1 <- ggsurvplot(km_fit, data = df,
pval = TRUE,
pval.coord = c(10, 0.1),
conf.int = FALSE,
risk.table = TRUE,
risk.table.height = 0.2,
risk.table.y.text = TRUE,
risk.table.title = "Number at Risk (number censored)",
xlim = c(0, 28),
break.time.by = 7,
legend.title = " ",
palette = c("#E63946", "#457B9D"),
ggtheme = theme_minimal(),
legend.labs = c("High", "Low")
)
Figure_6.1 <- Figure_6.1 + ggtitle("28-day Survival Curves for ADRB2")
# Combine the plots (ADR1 and ADR2) side by side with the risk tables
Figure_6_combined <- plot_grid(Figure_6$plot, Figure_6$table, ncol = 1, align = "v", rel_heights = c(3, 1))
Figure_6.1_combined <- plot_grid(Figure_6.1$plot, Figure_6.1$table, ncol = 1, align = "v", rel_heights = c(3, 1))
# Remove titles from individual plots
Figure_6_combined <- Figure_6_combined + theme(plot.title = element_blank())
Figure_6.1_combined <- Figure_6.1_combined + theme(plot.title = element_blank())
# Combine both survival plots (ADR1 and ADR2) together with a title
Figure6_tot <- plot_grid(Figure_6_combined, Figure_6.1_combined, ncol = 2, align = "h")
# Add a title to the combined figure
Figure6_tot <- Figure6_tot + plot_annotation(title = "28-day Survival Curves for ADRB expressed in Lymphocytes")
Figure6_tot
df_surv_m1 <- df %>%
filter(!is.na(ADR1m_D0_median), !is.na(diff_days_J28), !is.na(Outcome_censored_28))
# ADRB1 on Monocytes
km_fit_m1 <- survfit(Surv(diff_days_J28, Outcome_censored_28) ~ ADR1m_D0_median, data = df_surv_m1)
Fig_m1 <- ggsurvplot(km_fit_m1, data = df,
pval = TRUE,
pval.coord = c(10, 0.1),
conf.int = FALSE,
risk.table = TRUE,
risk.table.height = 0.2,
risk.table.y.text = TRUE,
risk.table.title = "Number at Risk (number censored)",
xlim = c(0, 28),
break.time.by = 7,
legend.title = " ",
palette = c("#E63946", "#457B9D"),
ggtheme = theme_minimal(),
legend.labs = c("High", "Low")
)
Fig_m1 <- Fig_m1 + ggtitle("28-day Survival Curves for ADRB1")
# ADRB2 on Monocytes
km_fit_m2 <- survfit(Surv(diff_days_J28, Outcome_censored_28) ~ ADR2m_D0_median, data = df)
Fig_m2 <- ggsurvplot(km_fit_m2, data = df,
pval = TRUE,
pval.coord = c(10, 0.1),
conf.int = FALSE,
risk.table = TRUE,
risk.table.height = 0.2,
risk.table.y.text = TRUE,
risk.table.title = "Number at Risk (number censored)",
xlim = c(0, 28),
break.time.by = 7,
legend.title = " ",
palette = c("#E63946", "#457B9D"),
ggtheme = theme_minimal(),
legend.labs = c("High", "Low")
)
Fig_m2 <- Fig_m2 + ggtitle("28-day Survival Curves for ADRB2")
# Combine plots
Fig_m1_combined <- plot_grid(Fig_m1$plot, Fig_m1$table, ncol = 1, align = "v", rel_heights = c(3, 1))
Fig_m2_combined <- plot_grid(Fig_m2$plot, Fig_m2$table, ncol = 1, align = "v", rel_heights = c(3, 1))
Fig_m1_combined <- Fig_m1_combined + theme(plot.title = element_blank())
Fig_m2_combined <- Fig_m2_combined + theme(plot.title = element_blank())
Fig_monocyte_survival <- Fig_m1_combined + Fig_m2_combined +
plot_annotation(title = "28-day Survival Curves for ADRB expressed in Monocytes")
Fig_monocyte_survival
To assess the prognostic significance of baseline β-adrenergic receptor expression, we stratified patients based on the median value of ADRB1 and ADRB2 measured at ECMO initiation (D0). Survival analysis was performed using Kaplan–Meier curves and the log-rank test to compare 28-day survival between patients with “High” vs “Low” expression. Patients were censored at the time of ECMO weaning, heart transplantation or LVAD implantation.
Although the differences did not reach statistical significance, a trend was observed indicating that patients with higher ADRB expression in both lymphocytes and monocytes at ECMO initiation had lower 28-day survival probability. This pattern suggests a potential association between adrenergic receptor activation and worse short-term outcomes.
plot_cytokine <- function(df, cytokine_name) {
# Colomn name
col_J0 <- paste0(cytokine_name, "_J0")
col_J3_J5 <- paste0(cytokine_name, "_J3_J5")
col_JS <- paste0(cytokine_name, "_JS")
# 1. dataset
data_plot <- df %>%
select(ID, Outcome, all_of(c(col_J0, col_J3_J5, col_JS))) %>%
mutate(
Outcomes = factor(case_when(
Outcome %in% c("Bridge to LVAD", "Bridge to Transplant") ~ "Bridge to transplant or LVAD",
TRUE ~ as.character(Outcome)
), levels = c("Death", "Bridge to transplant or LVAD", "ECMO Weaning"))
) %>%
pivot_longer(
cols = c(all_of(c(col_J0, col_J3_J5, col_JS))),
names_to = "cyto_time",
values_to = "value"
) %>%
mutate(
time = case_when(
cyto_time == col_J0 ~ 1,
cyto_time == col_J3_J5 ~ 2,
cyto_time == col_JS ~ 3,
TRUE ~ NA_real_
),
value_log = log10(value)
)
# 2. n
count_data <- data_plot %>%
filter(!is.na(value)) %>%
group_by(cyto_time, Outcomes) %>%
summarise(n = n(), .groups = "drop") %>%
filter(!is.na(Outcomes))
# 3. Models
mod1 <- lmerTest::lmer(
value_log ~ time * Outcomes + (1 + time | ID),
data = data_plot
)
mod_no_interaction <- lmerTest::lmer(
value_log ~ time + Outcomes + (1 + time | ID),
data = data_plot
)
p_anova <- anova(mod_no_interaction, ddf = "Kenward-Roger")
p_time <- format.pval(p_anova["time", "Pr(>F)"], digits = 3, eps = .001)
p_outcome <- format.pval(p_anova["Outcomes", "Pr(>F)"], digits = 3, eps = .001)
# 4. Plot
p <- ggplot(data_plot, aes(x = cyto_time, y = value_log, fill = Outcomes)) +
geom_boxplot(position = position_dodge(0.7), width = 0.6) +
scale_y_log10() +
labs(
title = paste("Association between", cytokine_name, "and outcomes"),
x = "Times of measurement",
y = paste("Log10 of", cytokine_name, "(pg/mL)")
) +
geom_text(data = count_data, aes(
x = cyto_time,
y = -min(data_plot$value_log, na.rm = TRUE) * 0.01,
label = paste0("n=", n),
group = Outcomes
), position = position_dodge(0.7), size = 4, vjust = 1) +
scale_x_discrete(labels = setNames(
c("implantation", "day 3 to 5", "explantation"),
c(col_J0, col_J3_J5, col_JS)
)) +
scale_fill_manual(values = c(
"ECMO Weaning" = "#4DAF4A",
"Bridge to transplant or LVAD" = "#377EB8",
"Death" = "#F8766D"
)) +
annotate(
"text",
x = 1,
y = max(data_plot$value_log, na.rm = TRUE) * 1.1,
label = paste0("p[time] == ", p_time, "*';'*~p[outcomes] == ", p_outcome),
parse = TRUE,
size = 4,
hjust = 0
)
print(p)
}
plot_cytokine(df, "IL6")
IL-6 levels were lower in patients successfully weaned from ECMO, particularly at baseline (day 0) and at explantation. A significant effect of time was observed (p = 0.02), while the association with clinical outcomes showed a trend without reaching statistical significance (p = 0.055).
plot_cytokine_linear <- function(df, cytokine_name = "IL-33") {
col_J0 <- paste0(cytokine_name, "_J0")
col_J3_J5 <- paste0(cytokine_name, "_J3_J5")
col_JS <- paste0(cytokine_name, "_JS")
data_plot <- df %>%
select(ID, Outcome, all_of(c(col_J0, col_J3_J5, col_JS))) %>%
mutate(
Outcomes = factor(case_when(
Outcome %in% c("Bridge to LVAD", "Bridge to Transplant") ~ "Bridge to transplant or LVAD",
TRUE ~ as.character(Outcome)
), levels = c("Death", "Bridge to transplant or LVAD", "ECMO Weaning"))
) %>%
pivot_longer(cols = all_of(c(col_J0, col_J3_J5, col_JS)),
names_to = "cyto_time", values_to = "value") %>%
filter(!is.na(value)) %>%
mutate(
time = case_when(
cyto_time == col_J0 ~ 1,
cyto_time == col_J3_J5 ~ 2,
cyto_time == col_JS ~ 3,
TRUE ~ NA_real_
)
)
count_data <- data_plot %>%
group_by(cyto_time, Outcomes) %>%
summarise(n = n(), .groups = "drop") %>%
filter(!is.na(Outcomes))
mod1 <- lmerTest::lmer(value ~ time * Outcomes + (1 + time | ID), data = data_plot)
mod_no_interaction <- lmerTest::lmer(value ~ time + Outcomes + (1 + time | ID), data = data_plot)
p_anova <- anova(mod_no_interaction, ddf = "Kenward-Roger")
p_time <- format.pval(p_anova["time", "Pr(>F)"], digits = 3, eps = .001)
p_outcome <- format.pval(p_anova["Outcomes", "Pr(>F)"], digits = 3, eps = .001)
ggplot(data_plot, aes(x = cyto_time, y = value, fill = Outcomes)) +
geom_boxplot(position = position_dodge(0.7), width = 0.6) +
labs(
title = "Association between IL-33 and outcomes",
x = "Times of measurement",
y = "IL-33 (pg/mL)"
) +
geom_text(data = count_data, aes(
x = cyto_time,
y = min(data_plot$value, na.rm = TRUE) * 0.9,
label = paste0("n=", n),
group = Outcomes
), position = position_dodge(0.7), size = 4, vjust = 1) +
scale_x_discrete(labels = setNames(
c("implantation", "day 3 to 5", "explantation"),
c(col_J0, col_J3_J5, col_JS)
)) +
scale_fill_manual(values = c(
"ECMO Weaning" = "#4DAF4A",
"Bridge to transplant or LVAD" = "#377EB8",
"Death" = "#F8766D"
)) +
annotate("text", x = 1, y = max(data_plot$value, na.rm = TRUE) * 1.1,
label = paste0("p[time] == ", p_time, "*';'*~p[outcomes] == ", p_outcome),
parse = TRUE, size = 4, hjust = 0)
}
plot_cytokine_linear(df, "IL-33")
Plasma concentrations of IL-33 were undetectable in the vast majority of patients, with only three samples showing quantifiable values above the detection threshold. Due to this distribution, IL-33 levels were visualized on a linear scale rather than log-transformed.
In addition to IL-33, we also measured IL-1β, IL-2 and IL-4; however, these cytokines were consistently undetectable across all samples, preventing any meaningful quantitative analysis. Similarly, plasma concentrations of IL-12p70, IFN-γ and Granzyme B were below the detection threshold in all patients. IL-17A was undetectable in the vast majority of samples, with only one patient showing a quantifiable value above the lower limit of detection. TNF-α was quantifiable in only two patients. Due to the pervasive non-detectability of these markers, no statistical comparison across timepoints or clinical outcomes was feasible.
plot_cytokine(df, "IL5")
Among the cytokines measured, IL-5 concentrations appeared lower in patients successfully weaned from ECMO, although differences across timepoints and outcomes did not reach statistical significance.
plot_cytokine(df, "IL8/CXCL8")
Also IL-8 concentrations appeared clear lower in patients successfully weaned from ECMO, although differences across timepoints and outcomes did not reach statistical significance.
plot_cytokine <- function(df, cytokine_name) {
# Column names for the three timepoints
col_J0 <- paste0(cytokine_name, "_J0")
col_J3_J5 <- paste0(cytokine_name, "_J3_J5")
col_JS <- paste0(cytokine_name, "_JS")
# 1. Prepare dataset
data_plot <- df %>%
select(ID, Outcome, all_of(c(col_J0, col_J3_J5, col_JS))) %>%
mutate(
Outcomes = factor(case_when(
Outcome %in% c("Bridge to LVAD", "Bridge to Transplant") ~ "Bridge to transplant or LVAD",
TRUE ~ as.character(Outcome)
), levels = c("Death", "Bridge to transplant or LVAD", "ECMO Weaning"))
) %>%
pivot_longer(
cols = c(all_of(c(col_J0, col_J3_J5, col_JS))),
names_to = "cyto_time",
values_to = "value"
) %>%
mutate(
time = case_when(
cyto_time == col_J0 ~ 1,
cyto_time == col_J3_J5 ~ 2,
cyto_time == col_JS ~ 3,
TRUE ~ NA_real_
),
value_log = log10(value)
)
# 2. Count number per group
count_data <- data_plot %>%
filter(!is.na(value)) %>%
group_by(cyto_time, Outcomes) %>%
summarise(n = n(), .groups = "drop") %>%
filter(!is.na(Outcomes))
# 3. Fit models
mod1 <- lmerTest::lmer(
value_log ~ time * Outcomes + (1 + time | ID),
data = data_plot
)
mod_no_interaction <- lmerTest::lmer(
value_log ~ time + Outcomes + (1 + time | ID),
data = data_plot
)
p_anova <- anova(mod_no_interaction, ddf = "Kenward-Roger")
p_time_raw <- format.pval(p_anova["time", "Pr(>F)"], digits = 3, eps = .001)
p_outcome_raw <- format.pval(p_anova["Outcomes", "Pr(>F)"], digits = 3, eps = .001)
# Safe formatting for p-values with '<'
p_time <- if (grepl("<", p_time_raw)) "\"<2e-16\"" else p_time_raw
p_outcome <- if (grepl("<", p_outcome_raw)) "\"<2e-16\"" else p_outcome_raw
# 4. Plot
p <- ggplot(data_plot, aes(x = cyto_time, y = value_log, fill = Outcomes)) +
geom_boxplot(position = position_dodge(0.7), width = 0.6) +
scale_y_log10() +
labs(
title = paste("Association between", cytokine_name, "and outcomes"),
x = "Times of measurement",
y = paste("Log10 of", cytokine_name, "(pg/mL)")
) +
geom_text(data = count_data, aes(
x = cyto_time,
y = -min(data_plot$value_log, na.rm = TRUE) * 0.01,
label = paste0("n=", n),
group = Outcomes
), position = position_dodge(0.7), size = 4, vjust = 1) +
scale_x_discrete(labels = setNames(
c("implantation", "day 3 to 5", "explantation"),
c(col_J0, col_J3_J5, col_JS)
)) +
scale_fill_manual(values = c(
"ECMO Weaning" = "#4DAF4A",
"Bridge to transplant or LVAD" = "#377EB8",
"Death" = "#F8766D"
)) +
annotate(
"text",
x = 1,
y = max(data_plot$value_log, na.rm = TRUE) * 1.1,
label = paste0("p[time] == ", p_time, "*';'*~p[outcomes] == ", p_outcome),
parse = TRUE,
size = 4,
hjust = 0
)
print(p)
}
plot_cytokine(df, "IL10")
IL-10 levels showed a statistically significant change over time (p < 0.05), indicating dynamic modulation during the clinical course. However, no significant differences were observed across the defined outcome groups.
plot_cytokine(df, "GDF-15")
GDF-15 concentrations appeared lower in patients successfully weaned from ECMO, although differences across timepoints and outcomes did not reach statistical significance.
Not all cytokines, whether pro-inflammatory or anti-inflammatory, were consistently measurable in our cohort. Based on the data, cytokine levels do not appear to be significantly associated with clinical outcomes. Interestingly, both pro-inflammatory and anti-inflammatory cytokines tend to decrease over time.
anti_infiamm <- c("IL10", "GDF-15")
pro_infiamm <- c("IL5", "IL6", "IL8/CXCL8")
plot_cytokines_group <- function(df, timepoint = "J3_J5") {
cols_anti <- paste0(anti_infiamm, "_", timepoint)
cols_pro <- paste0(pro_infiamm, "_", timepoint)
df_sub <- df %>%
select(ID, Outcome, all_of(c(cols_anti, cols_pro))) %>%
mutate(
Outcomes = factor(case_when(
Outcome %in% c("Bridge to LVAD", "Bridge to Transplant") ~ "Bridge to transplant or LVAD",
TRUE ~ Outcome
), levels = c("Death", "Bridge to transplant or LVAD", "ECMO Weaning"))
)
# Reshape long
df_long <- df_sub %>%
pivot_longer(cols = -c(ID, Outcome, Outcomes),
names_to = "cytokine_time", values_to = "value") %>%
mutate(
cytokine = sub(paste0("_", timepoint), "", cytokine_time),
# Negativ sign for anti-inflamm on left
value_plot = ifelse(cytokine %in% anti_infiamm, -value, value),
cytokine = factor(cytokine, levels = c(anti_infiamm, pro_infiamm))
) %>%
filter(!is.na(value))
# Colors
colori <- c(
"ECMO Weaning" = "#4DAF4A",
"Bridge to transplant or LVAD" = "#377EB8",
"Death" = "#F8766D"
)
# Plot
ggplot(df_long, aes(x = value_plot, y = fct_rev(cytokine), fill = Outcomes)) +
geom_bar(stat = "summary", fun = mean, position = "dodge", width = 0.7) +
scale_fill_manual(values = colori) +
scale_x_continuous(
breaks = scales::pretty_breaks(n = 10),
labels = abs,
name = "Mean cytokine level (pg/mL)"
) +
labs(title = paste("Cytokine levels at", timepoint, "by outcome group"),
y = "Cytokine") +
theme_minimal() +
theme(
legend.position = "right"
)
}
plot_cytokines_group(df, "J0")
plot_cytokines_group(df, "J3_J5")
plot_cytokines_group(df, "JS")
anti_infiamm <- c("IL10", "GDF-15")
pro_infiamm <- c("IL5", "IL6", "IL8/CXCL8")
plot_cytokines_group_log_vertical <- function(df, timepoint = "J3_J5") {
cols_anti <- paste0(anti_infiamm, "_", timepoint)
cols_pro <- paste0(pro_infiamm, "_", timepoint)
df_sub <- df %>%
select(ID, Outcome, all_of(c(cols_anti, cols_pro))) %>%
mutate(
Outcomes = factor(case_when(
Outcome %in% c("Bridge to LVAD", "Bridge to Transplant") ~ "Bridge to transplant or LVAD",
TRUE ~ Outcome
), levels = c("Death", "Bridge to transplant or LVAD", "ECMO Weaning"))
)
df_long <- df_sub %>%
pivot_longer(cols = -c(ID, Outcome, Outcomes),
names_to = "cytokine_time", values_to = "value") %>%
mutate(
cytokine = sub(paste0("_", timepoint), "", cytokine_time),
# Replace zeros or negatives before log10
value_adj = ifelse(value <= 0, 1e-7, value),
value_log = log10(value_adj),
# Group label
group = case_when(
cytokine %in% anti_infiamm ~ "Anti-inflammatory",
cytokine %in% pro_infiamm ~ "Pro-inflammatory",
TRUE ~ "Other"
),
# Custom sign flipping: IL10 positive, IL5 negative, other anti-inflam negative, pro positive
value_plot = case_when(
cytokine == "IL10" ~ value_log,
cytokine == "IL5" ~ -value_log,
group == "Anti-inflammatory" & cytokine != "IL10" ~ -value_log,
TRUE ~ value_log
),
cytokine = factor(cytokine, levels = c(anti_infiamm, pro_infiamm)),
group = factor(group, levels = c("Anti-inflammatory", "Pro-inflammatory"))
) %>%
filter(!is.na(value_log))
colori <- c(
"ECMO Weaning" = "#4DAF4A",
"Bridge to transplant or LVAD" = "#377EB8",
"Death" = "#F8766D"
)
# Plot base
p <- ggplot(df_long, aes(x = value_plot, y = fct_rev(cytokine), fill = Outcomes)) +
geom_bar(stat = "summary", fun = mean, position = position_dodge(width = 0.7), width = 0.6) +
scale_fill_manual(values = colori) +
scale_x_continuous(
breaks = scales::pretty_breaks(n = 10),
labels = function(x) round(abs(x), 2),
name = "Mean log10 cytokine level (pg/mL)"
) +
labs(title = paste("Log10 Cytokine levels at", timepoint, "by outcome group"),
y = "Cytokine") +
theme_minimal(base_size = 14) +
theme(
legend.position = "right",
axis.text.y = element_text(size = 12),
plot.margin = margin(t = 35, r = 20, b = 10, l = 20)
)
# Add top annotation for anti/pro groups with grid::annotation_custom
library(grid)
anti_label <- textGrob("Anti-inflammatory", x = unit(0.15, "npc"), y = unit(0.98, "npc"),
gp = gpar(fontsize = 14, fontface = "bold"))
pro_label <- textGrob("Pro-inflammatory", x = unit(0.85, "npc"), y = unit(0.98, "npc"),
gp = gpar(fontsize = 14, fontface = "bold"))
p + annotation_custom(anti_label) + annotation_custom(pro_label)
}
plot_cytokines_group_log_vertical(df, "J0")
At day 0, patients who later achieved ECMO weaning tended to have higher levels of anti-inflammatory cytokines such as IL-10. GDF-15 levels were comparable between patients who died and those who were weaned. In contrast, pro-inflammatory cytokines—particularly IL-6 and IL-8—were elevated in patients who died.
(For Antoine: For better graphical representation, I choose a logarithmic scale. Among all timepoints, I will focused primarily on day 0, aiming to identify potential predictors or markers that could help clinicians recognize patients with a more aggressive disease course and possibly guide management decisions, although currently no targeted therapy has been proven effective.)
# CK dosable
cyto_selected <- c("IL10", "GDF-15", "IL5", "IL6", "IL8/CXCL8")
anti_infiamm <- c("IL10", "GDF-15")
# At J0
cyto_cols <- paste0(cyto_selected, "_J0")
# Cutoff mediana + indosable
quantile_3groups_label_simple <- function(x, detection_limit = 0.001) {
x_label <- ifelse(is.na(x) | x < detection_limit, "Indosabile", NA_character_)
x_valid <- x[!(is.na(x) | x < detection_limit)]
med_val <- median(x_valid, na.rm = TRUE)
x_label[!(is.na(x) | x < detection_limit)] <- ifelse(x_valid <= med_val, "Low", "High")
factor(x_label, levels = c("Indosabile", "Low", "High"))
}
# Groups High/Low ADRB1 e ADRB2 in lympho e mono (median)
df <- df %>%
mutate(
ADR1m_group = ifelse(ADR1m_J0 > median(ADR1m_J0, na.rm = TRUE), "High", "Low"),
ADR2m_group = ifelse(ADR2m_J0 > median(ADR2m_J0, na.rm = TRUE), "High", "Low"),
ADR1_group = ifelse(ADR1_J0 > median(ADR1_J0, na.rm = TRUE), "High", "Low"),
ADR2_group = ifelse(ADR2_J0 > median(ADR2_J0, na.rm = TRUE), "High", "Low")
)
df_quart <- df %>%
select(ID, all_of(cyto_cols), ADR1m_group, ADR2m_group, ADR1_group, ADR2_group) %>%
mutate(across(all_of(cyto_cols), ~ quantile_3groups_label_simple(.))) %>%
filter(!is.na(ADR1m_group), !is.na(ADR2m_group), !is.na(ADR1_group), !is.na(ADR2_group))
make_long <- function(df, group_col, cell_type, adr_type) {
df %>%
select(ID, all_of(cyto_cols), all_of(group_col)) %>%
pivot_longer(cols = all_of(cyto_cols), names_to = "cytokine", values_to = "level") %>%
mutate(
cytokine = sub("_J0$", "", cytokine),
adr_group = .data[[group_col]],
cell_type = cell_type,
adr_type = adr_type
) %>%
select(-all_of(group_col))
}
df_long_lymph_ADR1 <- make_long(df_quart, "ADR1_group", "Lymphocytes", "ADRB1")
df_long_lymph_ADR2 <- make_long(df_quart, "ADR2_group", "Lymphocytes", "ADRB2")
df_long_mono_ADR1 <- make_long(df_quart, "ADR1m_group", "Monocytes", "ADRB1")
df_long_mono_ADR2 <- make_long(df_quart, "ADR2m_group", "Monocytes", "ADRB2")
heatmap_data <- bind_rows(df_long_lymph_ADR1, df_long_lymph_ADR2, df_long_mono_ADR1, df_long_mono_ADR2) %>%
mutate(
cytokine = factor(cytokine, levels = cyto_selected),
adr_group = factor(adr_group, levels = c("Low", "High")),
cell_type = factor(cell_type, levels = c("Lymphocytes", "Monocytes")),
adr_type = factor(adr_type, levels = c("ADRB1", "ADRB2"))
)
palette_colors <- c(
"Indosabile" = "white",
"Low" = "forestgreen",
"High" = "firebrick"
)
# Plot
ggplot(heatmap_data, aes(x = adr_group, y = cytokine, fill = quartile)) +
# Heatmap dati
geom_tile(data = heatmap_data, aes(x = adr_group, y = cytokine, fill = level), color = "grey80", width = 0.8) +
facet_grid(cell_type ~ adr_type) +
scale_fill_manual(values = palette_colors, na.value = "grey90") +
scale_x_discrete(expand = expansion(add = c(1, 0.1))) +
coord_cartesian(clip = "off") +
labs(
title = "Heatmap of Selected Cytokine Levels at Day 0 by ADRB Expression Groups",
x = "ADRB Expression Group",
y = "Cytokine",
fill = "Level"
) +
theme_minimal() +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(color = "black", angle = 45, hjust = 1),
panel.grid = element_blank(),
strip.background = element_blank(),
strip.text.x = element_text(hjust = 0.5),
strip.text.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5)
)
# I don't understand why IL5 and IL8 are white